Supporting code, figures and tables for
Urban birds’ flight responses were largely unaffected by the COVID-19 shutdowns
Peter Mikula, Martin Bulla, Daniel T. Blumstein, Yanina Benedetti, Kristina Floigl, Jukka Jokimäki, Marja-Liisa Kaisanlahti-Jokimäki, Gábor Markó, Federico Morelli, Anders Pape Møller, Anastasiia Siretckaia, Sára Szakony, Michael A. Weston, Farah Abou Zeid, Piotr Tryjanowski & Tomáš Albrecht
created by Martin Bulla
2023-06-26 14:49:47
Code to load tools and data
# packages
require(arm)
require(data.table)
require(effects)
require(foreach)
require(ggimage)
require(ggplot2)
require(ggpubr)
require(ggsci)
require(ggtext)
require(grid)
require(gtable)
require(here)
require(kableExtra)
require(MASS)
require(multcomp)
require(optimx)
require(performance)
require(PerformanceAnalytics)
require(png)
require(RColorBrewer)
require(rmeta)
require(rphylopic)
require(scales)
require(viridis)
# constants
save_plot = TRUE
round_ = 3 # number of decimal places to round model coefficients
nsim = 5000 # number of simulations to extract estimates and 95%CrI
ax_lines = "grey60" # defines color of the axis lines
#colors <- c("#999999", "#E69F00", "#56B4E9") #viridis(3)
set.seed(42)
#width_ = .7 # spacing between error bars
#col_ = c(brewer.pal(n =12, name = "Paired"), 'grey30','grey80')
# functions
# to add images to panels
annotation_custom2 <- function (grob, xmin = -Inf, xmax = Inf, ymin = -Inf, ymax = Inf, data) {
layer(data = data, stat = StatIdentity, position = PositionIdentity,
geom = ggplot2:::GeomCustomAnn,
inherit.aes = TRUE, params = list(grob = grob,
xmin = xmin, xmax = xmax,
ymin = ymin, ymax = ymax))
}
# to remove ggplot components
gtable_filter_remove <- function (x, name, trim = TRUE){
matches <- !(x$layout$name %in% name)
x$layout <- x$layout[matches, , drop = FALSE]
x$grobs <- x$grobs[matches]
if (trim)
x <- gtable_trim(x)
x
}
# customized ggplot theme
theme_MB = theme(
title = element_text(size=8, colour="grey30"),
axis.line = element_blank(),
#axis.line = element_line(colour="grey70", size=0.25),
axis.title = element_text(size=7, colour="grey30"),
axis.title.y = element_text(vjust=3.5),
axis.title.x = element_text(vjust=1),
axis.text = element_text(size=6),#, vjust = 0.5, hjust=1),# margin=units(0.5,"mm")),
axis.ticks.length=unit(0.5,"mm"),
axis.ticks = element_line(colour = "grey70", size = 0.1),
#axis.ticks.margin,
strip.text.x = element_text(size = 6, color="grey30", margin=margin(1,1,1,1,"mm")), #grey50
strip.text.y = element_text(size = 6, color="grey30", margin=margin(1,1,1,1,"mm")), #grey50
strip.background = element_rect(fill="grey99",colour="grey70", size=0.25),
#strip.background = element_blank(),
#strip.text = element_blank(),
panel.spacing = unit(0, "mm"),
panel.background=element_blank(),
panel.border = element_rect(colour="grey70", size=0.1, fill = NA), #panel.border=element_blank(),
panel.grid = element_blank(),
legend.text=element_text(size=6),
legend.title=element_text(size=6),
legend.key = element_rect(colour = NA, fill = NA),
legend.key.height= unit(0.5,"line"),
legend.key.width = unit(0.25, "cm"),
legend.margin = margin(0,0,0,0, unit="cm"),
legend.box.margin = margin(l = -6), #legend.justification = c(-1,0),
legend.background = element_blank()
)
# for estimates
est_out =function(model = m, label = "", nsim = 5000){
bsim = sim(model, n.sim=5000)
v = apply(bsim@fixef, 2, quantile, prob=c(0.5))
ci = apply(bsim@fixef, 2, quantile, prob=c(0.025,0.975))
sd = apply(bsim@fixef, 2, sd)
data.table(predictor=rownames(coef(summary(model))),estimate=v, lwr=ci[1,], upr=ci[2,], sd = sd, model = paste(label, "N =", nobs(model)))
}
# change color
change_col = function(replace_black, theimg) {
r_b = col2rgb(replace_black) / 255
#theimg[theimg == 1] <- 2
for (i in 1:3) {
theimg[,,i][theimg[,,i] == 0] <- r_b[i]
}
return(theimg)
}
# for Supplementary Table output based on sim
m_out = function(model = m, type = "mixed",
name = "define", dep = "define", fam = 'Gaussian',
round_ = 3, nsim = 5000, aic = FALSE, save_sim = here::here('Data/model_sim/'), back_tran = FALSE, perc_ = 1){
# perc_ 1 = proportion or 100%
bsim = sim(model, n.sim=nsim)
if(save_sim!=FALSE){save(bsim, file = paste0(save_sim, name,'.RData'))}
if(type != "mixed"){
v = apply(bsim@coef, 2, quantile, prob=c(0.5))
ci = apply(bsim@coef, 2, quantile, prob=c(0.025,0.975))
if(back_tran == TRUE & fam == "binomial"){
v = perc_*plogis(v)
ci = perc_*plogis(ci)
}
if(back_tran == TRUE & fam == "binomial_logExp"){
v = perc_*(1-plogis(v))
ci = perc_*(1-plogis(ci))
ci = rbind(ci[2,],ci[1,])
}
if(back_tran == TRUE & fam == "Poisson"){
v = exp(v)
ci = exp(ci)
}
oi=data.frame(type='fixed',effect=rownames(coef(summary(model))),estimate=v, lwr=ci[1,], upr=ci[2,])
rownames(oi) = NULL
oi$estimate_r=round(oi$estimate,round_)
oi$lwr_r=round(oi$lwr,round_)
oi$upr_r=round(oi$upr,round_)
if(perc_ == 100){
oi$estimate_r = paste0(oi$estimate_r,"%")
oi$lwr_r = paste0(oi$lwr_r,"%")
oi$upr_r = paste0(oi$upr_r,"%")
}
x=data.table(oi[c('type',"effect", "estimate_r","lwr_r",'upr_r')])
}else{
v = apply(bsim@fixef, 2, quantile, prob=c(0.5))
ci = apply(bsim@fixef, 2, quantile, prob=c(0.025,0.975))
if(back_tran == TRUE & fam == "binomial"){
v = perc_*plogis(v)
ci = perc_*plogis(ci)
}
if(back_tran == TRUE & fam == "binomial_logExp"){
v = perc_*(1-plogis(v))
ci = perc_*(1-plogis(ci))
ci = rbind(ci[2,],ci[1,])
}
if(back_tran == TRUE & fam == "Poisson"){
v = exp(v)
ci = exp(ci)
}
oi=data.table(type='fixed',effect=rownames(coef(summary(model))),estimate=v, lwr=ci[1,], upr=ci[2,])
rownames(oi) = NULL
oi[,estimate_r := round(estimate,round_)]
oi[,lwr_r := round(lwr,round_)]
oi[,upr_r :=round(upr,round_)]
if(perc_ == 100){
oi[,estimate_r := paste0(estimate_r,"%")]
oi[,lwr_r := paste0(lwr_r,"%")]
oi[,upr_r := paste0(upr_r,"%")]
}
oii=oi[,c('type',"effect", "estimate_r","lwr_r",'upr_r')]
l=data.frame(summary(model)$varcor)
l = l[is.na(l$var2),]
l$var1 = ifelse(is.na(l$var1),"",l$var1)
l$pred = paste(l$grp,l$var1)
q050={}
q025={}
q975={}
pred={}
# variance of random effects
for (ran in names(bsim@ranef)) {
ran_type = l$var1[l$grp == ran]
for(i in ran_type){
q050=c(q050,quantile(apply(bsim@ranef[[ran]][,,ran_type], 1, var), prob=c(0.5)))
q025=c(q025,quantile(apply(bsim@ranef[[ran]][,,ran_type], 1, var), prob=c(0.025)))
q975=c(q975,quantile(apply(bsim@ranef[[ran]][,,ran_type], 1, var), prob=c(0.975)))
pred= c(pred,paste(ran, i))
}
}
# residual variance
q050=c(q050,quantile(bsim@sigma^2, prob=c(0.5)))
q025=c(q025,quantile(bsim@sigma^2, prob=c(0.025)))
q975=c(q975,quantile(bsim@sigma^2, prob=c(0.975)))
pred= c(pred,'Residual')
ri=data.table(type='random',effect=pred, estimate_r=round(100*q050/sum(q050)), lwr_r=round(100*q025/sum(q025)), upr_r=round(100*q975/sum(q975)))
ri[lwr_r>upr_r, lwr_rt := upr_r]
ri[lwr_r>upr_r, upr_rt := lwr_r]
ri[!is.na(lwr_rt), lwr_r := lwr_rt]
ri[!is.na(upr_rt), upr_r := upr_rt]
ri$lwr_rt = ri$upr_rt = NULL
ri[,estimate_r := paste0(estimate_r,'%')]
ri[,lwr_r := paste0(lwr_r,'%')]
ri[,upr_r := paste0(upr_r,'%')]
x = data.table(rbind(oii,ri))
}
x[1, model := name]
x[1, response := dep]
x[1, error_structure := fam]
N = length(resid(model))
x[1, N := N ]
x=x[ , c('model', 'response', 'error_structure', 'N', 'type',"effect", "estimate_r","lwr_r",'upr_r')]
if (aic == TRUE){
x[1, AIC := AIC(update(model,REML = FALSE))]
}
if (aic == "AICc"){
aicc = AICc(model)
x[1, AICc := aicc]
}
if(type == "mixed" & nrow(x[type=='random' & estimate_r =='0%'])==0){
R2_mar = as.numeric(r2_nakagawa(model)$R2_marginal)
R2_con = as.numeric(r2_nakagawa(model)$R2_conditional)
x[1, R2_mar := R2_mar]
x[1, R2_con := R2_con]
}
x[is.na(x)] = ""
return(x)
}
# model assumption function
m_ass = function(name = 'define', mo = m0, dat = d, fixed = NULL, categ = NULL, trans = "none", spatial = TRUE, temporal = TRUE, PNG = TRUE, outdir = 'outdir', n_col=8, width_ = 10, height_ = 5.5){
l=data.frame(summary(mo)$varcor)
l = l[is.na(l$var2),]
nt = if(temporal==TRUE){1}else{0}
ns = if(spatial==TRUE){7}else{0}
n = 3+nrow(l)-1+length(fixed)+length(categ) + nt + ns
if(PNG == TRUE){
png(paste(outdir,name, ".png", sep=""), width=width_,height=height_,units="in",res=600) # width = 6
par(mfrow=c(4, n_col),tcl = -0.08, cex = 0.5, cex.main = 0.9,#ceiling(n/n_col),n_col)
oma = c(1,1,2,1),
mar = c(2, 2, 2, 1), mgp=c(1,0,0)
)
}else{
dev.new(width=width_,height=height_)
par(mfrow=c(4,n_col), tcl = -0.08, cex = 0.5, cex.main = 0.9,#ceiling(n/n_col),n_col)
oma = c(1,1,2,1),
mar = c(2, 2, 2, 1), mgp=c(1,0,0)
)
}
scatter.smooth(fitted(mo),resid(mo),col='grey');abline(h=0, lty=2, col ='red')
scatter.smooth(fitted(mo),sqrt(abs(resid(mo))), col='grey')
qqnorm(resid(mo), main=list("Normal Q-Q Plot: residuals"),col='grey');qqline(resid(mo), col = 'red')
#unique(l$grp[l$grp!="Residual"])
for(i in unique(l$grp[l$grp!="Residual"])){
#i = "mean_year"
ll=ranef(mo)[names(ranef(mo))==i][[1]]
if(ncol(ll)==1){
qqnorm(ll[,1], main = paste(i,names(ll)[1]),col='grey',);qqline(ll[,1], col ='red')
}else{
qqnorm(ll[,1], main = paste(i,names(ll)[1]),col='grey');qqline(ll[,1], col ='red')
qqnorm(ll[,2], main = paste(i,names(ll)[2]),col='grey');qqline(ll[,2], col ='red')
}
}
# variables
scatter={}
for (i in rownames(summary(mo)$coef)) {
#i = "lat_abs"
j=sub("\\).*", "", sub(".*\\(", "",i))
scatter[length(scatter)+1]=j
}
x = data.frame(scatter=unique(scatter)[2:length(unique(scatter))],
log_ = grepl("log",rownames(summary(mo)$coef)[2:length(unique(scatter))]), stringsAsFactors = FALSE)
for (i in 1:length(fixed)){
jj =fixed[i]
variable=dat[, ..jj][[1]]
if(trans[i]=='log'){
scatter.smooth(resid(mo)~log(variable),xlab=paste('log(',jj,')',sep=''), col = 'grey');abline(h=0, lwd=1, lty = 2, col ='red')
}else if(trans[i]=='abs'){
scatter.smooth(resid(mo)~abs(variable),xlab=paste('abs(',jj,')',sep=''), col = 'grey');abline(h=0, lwd=1, lty = 2, col ='red')
}else if(trans[i]=='sin'){scatter.smooth(resid(mo)~sin(variable),xlab=paste('sin(',jj,')',sep=''), col = 'grey');abline(h=0, lwd=1, lty = 2, col ='red')
}else if(trans[i]=='cos'){scatter.smooth(resid(mo)~cos(variable),xlab=paste('cos(',jj,')',sep=''), col = 'grey');abline(h=0, lwd=1, lty = 2, col ='red')
}else{
scatter.smooth(resid(mo)~variable,xlab=jj,col = 'grey');abline(h=0, lwd=1, lty = 2, col ='red')
}
}
if(length(categ)>0){
for(i in categ){
variable=dat[, ..i][[1]]
boxplot(resid(mo)~variable, medcol='grey', whiskcol='grey', staplecol='grey', boxcol='grey', outcol='grey', xlab = i);abline(h=0, lty=3, lwd=1, col = 'red')
}
}
if(temporal == TRUE){
acf(resid(mo), type="p", main=list("Temporal autocorrelation:\npartial series residual"))
}
if(spatial == TRUE){
spdata=data.frame(resid=resid(mo), x=dat$Lon, y=dat$Lat)
spdata$col=ifelse(spdata$resid<0,rgb(83,95,124,100, maxColorValue = 255),ifelse(spdata$resid>0,rgb(253,184,19,100, maxColorValue = 255), 'red'))
#cex_=c(1,2,3,3.5,4)
cex_=c(1,1.5,2,2.5,3)
spdata$cex=as.character(cut(abs(spdata$resid), 5, labels=cex_))
plot(spdata$x, spdata$y,col=spdata$col, cex=as.numeric(spdata$cex), pch= 16, main=list('Spatial distribution of residuals', cex=0.8), xlab = 'Longitude', ylab = 'Latitude')
legend("topright", pch=16, legend=c('>0','<0'), ,col=c(rgb(83,95,124,100, maxColorValue = 255),rgb(253,184,19,100, maxColorValue = 255)))
plot(spdata$x[spdata$resid<0], spdata$y[spdata$resid<0],col=spdata$col[spdata$resid<0], cex=as.numeric(spdata$cex[spdata$resid<0]), pch= 16, main=list('residuals <0'), xlab = 'Longitude', ylab = 'Latitude')
plot(spdata$x[spdata$resid>=0], spdata$y[spdata$resid>=0],col=spdata$col[spdata$resid>=0], cex=as.numeric(spdata$cex[spdata$resid>=0]), pch= 16, main=list('residual >=0'), xlab = 'Longitude', ylab = 'Latitude')
# EU
dat$res = resid(mo)
spdata=data.frame(resid = dat$res[dat$Country!='Australia'], x=dat$Lon[dat$Country!='Australia'], y=dat$Lat[dat$Country!='Australia'])
spdata$col=ifelse(spdata$resid<0,rgb(83,95,124,100, maxColorValue = 255),ifelse(spdata$resid>0,rgb(253,184,19,100, maxColorValue = 255), 'red'))
#cex_=c(1,2,3,3.5,4)
cex_=c(1,1.5,2,2.5,3)
spdata$cex=as.character(cut(abs(spdata$resid), 5, labels=cex_))
plot(spdata$x[spdata$resid<0], spdata$y[spdata$resid<0],col=spdata$col[spdata$resid<0], cex=as.numeric(spdata$cex[spdata$resid<0]), pch= 16, main=list('EU - residuals <0'), xlab = 'Longitude', ylab = 'Latitude')
plot(spdata$x[spdata$resid>=0], spdata$y[spdata$resid>=0],col=spdata$col[spdata$resid>=0], cex=as.numeric(spdata$cex[spdata$resid>=0]), pch= 16, main=list('EU residuals >=0)'), xlab = 'Longitude', ylab = 'Latitude')
# Australia
spdata=data.frame(resid = dat$res[dat$Country=='Australia'], x=dat$Lon[dat$Country=='Australia'], y=dat$Lat[dat$Country=='Australia'])
spdata$col=ifelse(spdata$resid<0,rgb(83,95,124,100, maxColorValue = 255),ifelse(spdata$resid>0,rgb(253,184,19,100, maxColorValue = 255), 'red'))
#cex_=c(1,2,3,3.5,4)
cex_=c(1,1.5,2,2.5,3)
spdata$cex=as.character(cut(abs(spdata$resid), 5, labels=cex_))
plot(spdata$x[spdata$resid<0], spdata$y[spdata$resid<0],col=spdata$col[spdata$resid<0], cex=as.numeric(spdata$cex[spdata$resid<0]), pch= 16, main=list('Australia residuals <0'), xlab = 'Longitude', ylab = 'Latitude')
plot(spdata$x[spdata$resid>=0], spdata$y[spdata$resid>=0],col=spdata$col[spdata$resid>=0], cex=as.numeric(spdata$cex[spdata$resid>=0]), pch= 16, main=list('Australia residuals >=0'), xlab = 'Longitude', ylab = 'Latitude')
}
mtext(stringr::str_wrap(paste(paste0(name," model: "), slot(mo,"call")[1],'(',slot(mo,"call")[2],sep=''), width = ceiling(nchar(paste(slot(mo,"call")[1],'(',slot(mo,"call")[2],sep=''))/2)+10), side = 3, line = 0, cex=0.5,outer = TRUE, col = 'darkblue') #ceiling(nchar(paste(slot(mo,"call")[1],'(',slot(mo,"call")[2],sep=''))/2)
if(PNG==TRUE){dev.off()}
}
# data
ph = fread(here::here('Data/phylopic.txt'))
setnames(ph, old = c('Name', 'Code'), new = c('genus2', 'uid'))
t = fread(here::here('Data/taxonomy.txt'))
g = fread(here::here('Data/google_mobility.txt')) #fwrite(d, here::here('Data/data.txt'), sep ='\t')
g[, Year := as.integer(substring(date, nchar(date)-3, nchar(date)))]
g[nchar(date)==9, date:=paste0('0',date)]
g[, date_ :=as.Date(date, format = '%d.%m.%Y')]
g[, Day :=yday(date_)]
g[country_region!='Australia', Day := Day-92 +1] # 1 April = start of breeding season (1st day) = 92 day of the year
g[country_region=='Australia', Day := Day-228 +1] # 15 Augusst = start of breeding season (1st day) = 228 day of the year
setnames(g, old = 'country_region', new ='Country')
d = fread(here::here('Data/data_corrected.txt')) #fwrite(d, here::here('Data/data.txt'), sep ='\t')
# add data and weekdays
x = fread(here::here("Data/date.txt"))[, .(IDObs, Date_corr)]
x[, date_:=as.Date(Date_corr, "%d.%m.%Y" )]
x[, weekday := weekdays(date_)]
d = merge(d, x[, .(IDObs, date_, weekday)], by = "IDObs")
#d[, date_ := as.Date(Day, origin = paste0(Year, '-01-01'))]
# adjust correct assignment of season (Year) for Australia
d[Country == 'Australia' & Year == 2020 & Covid == 0, Year:=2019]
d[Country == 'Australia' & Year == 2021 & Day>139, Year:=2020]
d[Country == 'Australia' & Year == 2022 & Day>139, Year:=2021]
d = d[order(Year, IDLocality, Day, Hour)]
d[Country %in% c("Czech_Republic", "Czech Republic"), Country := "Czechia"]
d[, genus := sub("_.*", "", Species)]
d[, sp_day_year := paste(Year, Species, Day, sep="_")]
d[, sp_loc := paste(Species, IDLocality, sep="_")]
d[, sp_country := paste(Species, Country, sep="\n")]
d[, rad:=(2*pi*Hour) / 24]
d[, Day_:= Day]
#d[, FID_z := scale(FID), by = Species]
#d[, SD_z := scale(SD), by = Species]
d[, FID_ln := log(FID)]
d[, SD_ln := log(SD)]
d[, body_ln := log(BodyMass)]
d[, flock_ln := log(FlockSize)]
d[, weekday := weekdays(date_)]
#d[Country == 'Australia', Day_:= abs(Day - 189)]
d1 = d[Covid == 1, .N, by = Species]
d2 = d[Covid == 0, .N, by = Species]
setnames(d1, old = 'N', new ='N_during')
setnames(d2, old = 'N', new ='N_before')
dd = merge(d1,d2) # species with data before and during
da = merge(d1,d2, all = TRUE)
d1p = d[Country!='Poland' & Covid == 1, .N, by = Species]
d2p = d[Country!='Poland' & Covid == 0, .N, by = Species]
setnames(d1p, old = 'N', new ='N_during')
setnames(d2p, old = 'N', new ='N_before')
ddp = merge(d1p,d2p) # species with data before and during, but without Poland data
p1 = d[Covid == 1, .N, by = .(IDLocality, Species)]
p2 = d[Covid == 0, .N, by = .(IDLocality, Species)]
setnames(p1, old = 'N', new ='N_during')
setnames(p2, old = 'N', new ='N_before')
pp = merge(p1,p2) # species-localities with data before and during
pa = merge(p1,p2, all = TRUE)
p1p = d[Country!='Poland' & Covid == 1, .N, by = .(IDLocality, Species)]
p2p = d[Country!='Poland' & Covid == 0, .N, by = .(IDLocality, Species)]
setnames(p1p, old = 'N', new ='N_during')
setnames(p2p, old = 'N', new ='N_before')
ppp = merge(p1p,p2p) # species-localities with data before and during,but without Poland data
p1p4 = d[Year!=2014 & Country!='Poland' & Covid == 1, .N, by = .(IDLocality, Species)]
p2p4 = d[Year!=2014 &Country!='Poland' & Covid == 0, .N, by = .(IDLocality, Species)]
setnames(p1p4, old = 'N', new ='N_during')
setnames(p2p4, old = 'N', new ='N_before')
ppp4 = merge(p1p4,p2p4) # species-localities with data before and during, but without Poland & 2014 data
s = d[Covid == 1]
s[, Nsp := .N, by ='Species']
s[, sp := gsub('[_]', ' ', Species)]
# add google mobility
s = merge(s, g[,.(Country, date_, parks_percent_change_from_baseline)], all.x = TRUE)
s[, year_ := as.character(Year)]
ss = s[!is.na(parks_percent_change_from_baseline)]
ss[, country_year := paste(Country, Year)] #table(paste(s$Country, s$Year))
ss[parks_percent_change_from_baseline<0, google := 'before_zero']
ss[parks_percent_change_from_baseline>0, google := 'after_zero']
ss[, sp_country_google:= paste(sp_country, google)]
g[, weekday := weekdays(date_)]Introduction
To facilitate transparency, the following document contains example code used to generate the results of the manuscript. Thus, apart from the supplementary figures and tables, below we display also main text figures. The figures and tables are ordered according to their appearance in the main text. The code is displayed upon clicking the “code” icon at the top right above each display item. Continuous variables were standardised by subtracting the mean and dividing by the standard deviation. For descriptions of variables see Methods of the paper. ‘Residual’ in tables indicates residual variance.
Effect sizes for Period
# predictions
# full
ms <- lmer(scale(log(FID)) ~
scale(log(SD)) +
scale(log(FlockSize)) +
scale(log(BodyMass)) +
scale(sin(rad)) + scale(cos(rad)) +
# scale(Day)+
scale(Temp) +
scale(Covid) +
(1 | Year) + (1 | weekday) + (1| genus) + (1| Species) + (1 | sp_day_year) + (scale(Covid) | Country) + (scale(Covid) | IDLocality) + (1 | sp_loc),
data = d, REML = FALSE,
control <- lmerControl(
optimizer = "optimx", optCtrl = list(method = "nlminb")
)
)
est_ms <- est_out(ms, "All: (1|Year) + (1|genus) + (1|Species) + (1|sp_day_year) + (scale(Covid)|Country) + (scale(Covid)|IDLocality) + (1|sp_loc)")
est_ms[, control_for_starting_distance := "yes"]
mx <- lmer(scale(log(FID)) ~
#scale(log(SD)) +
scale(log(FlockSize)) +
scale(log(BodyMass)) +
scale(sin(rad)) + scale(cos(rad)) +
# scale(Day)+
scale(Temp) +
scale(Covid) +
(1 | Year) + (1 | weekday) + (1 | genus) + (1 | Species) + (1 | sp_day_year) + (scale(Covid) | Country) + (scale(Covid) | IDLocality) + (1 | sp_loc),
data = d, REML = FALSE,
control = lmerControl(
optimizer ='optimx', optCtrl=list(method='nlminb'))
)
# (Covid|IDLocality) +
est_mx <- est_out(mx, "All: (1|Year) + (1|genus) + (1|Species) + (1|sp_day_year) + (scale(Covid)|Country) + (scale(Covid)|IDLocality) + (1|sp_loc)")
est_mx[, control_for_starting_distance := "no"]
# CZ - singular fits only due to genera estimated as zero (removing it changes no results)
cs <- lmer(scale(log(FID)) ~
scale(log(SD)) +
scale(log(FlockSize)) +
scale(log(BodyMass)) +
scale(sin(rad)) + scale(cos(rad)) +
# scale(Day)+
scale(Temp) +
scale(Covid) +
(1 | Year) + (1 | weekday) + (1 | genus) + (1| Species) + (1 | sp_day_year) + (scale(Covid)| IDLocality) + (1|sp_loc),
#(1 | Year) + (1 | weekday) + (1|genus) + (1 | Species) + (1 | sp_day_year) + (1 | IDLocality),
data = d[Country == "Czechia"], REML = FALSE
)
est_cs <- est_out(cs, "Czechia: (1|weekday)+(1|genus)+(1|Species)+(1|sp_day_year)+(scale(Covid)|IDLocality)+(1|sp_loc)")
est_cs[, control_for_starting_distance := "yes"]
cx <- lmer(scale(log(FID)) ~
#scale(log(SD)) +
scale(log(FlockSize)) +
scale(log(BodyMass)) +
scale(sin(rad)) + scale(cos(rad)) +
# scale(Day)+
scale(Temp) +
scale(Covid) +
(1 | Year) + (1 | weekday) + (1 | genus) + (1| Species) + (1 | sp_day_year) + (scale(Covid)| IDLocality) + (1|sp_loc),
#(1 | Year) + (1 | weekday) + (1 | genus) + (1 | Species) + (1 | sp_day_year) + (1 | IDLocality),
data = d[Country == "Czechia"], REML = FALSE
)
est_cx <- est_out(cx, "Czechia: (1|weekday)+(1|genus)+(1|Species)+(1|sp_day_year)+(scale(Covid)|IDLocality)+(1|sp_loc)")
est_cx[, control_for_starting_distance := "no"]
# FI
fs <- lmer(scale(log(FID)) ~
scale(log(SD)) +
scale(log(FlockSize)) +
scale(log(BodyMass)) +
scale(sin(rad)) + scale(cos(rad)) +
# scale(Day)+
scale(Temp) +
scale(Covid) +
(1 | Year) + (1 | weekday) + (1 | genus) + (1| Species) + (1 | sp_day_year) + (scale(Covid)| IDLocality) + (1|sp_loc),
data = d[Country == "Finland"], REML = FALSE
)
est_fs <- est_out(fs, "Finland: (1|weekday)+(1|genus)+(1|Species)+(1|sp_day_year)+(scale(Covid)|IDLocality)+(1|sp_loc)")
est_fs[, control_for_starting_distance := "yes"]
fx <- lmer(scale(log(FID)) ~
# scale(log(SD)) +
scale(log(FlockSize)) +
scale(log(BodyMass)) +
scale(sin(rad)) + scale(cos(rad)) +
# scale(Day)+
scale(Temp) +
scale(Covid) +
(1 | Year) + (1 | weekday) + (1 | genus) + (1| Species) + (1 | sp_day_year) + (scale(Covid)| IDLocality) + (1|sp_loc),
data = d[Country == "Finland"], REML = FALSE
)
est_fx <- est_out(fx, "Finland: (1|weekday)+(1|genus)+(1|Species)+(1|sp_day_year)+(scale(Covid)|IDLocality)+(1|sp_loc)")
est_fx[, control_for_starting_distance := "no"]
# HU - singular fits only due to sp_loc estimated as zero (removing it changes no results)
hs <- lmer(scale(log(FID)) ~
scale(log(SD)) +
scale(log(FlockSize)) +
scale(log(BodyMass)) +
scale(sin(rad)) + scale(cos(rad)) +
# scale(Day)+
scale(Temp) +
scale(Covid) +
(1 | Year) + (1 | weekday) + (1 | genus) + (1 | Species) + (1 | sp_day_year) + (scale(Covid) | IDLocality) + (1 | sp_loc),
# (1 | Year) + (1 | weekday) + (1|genus) + (1 | Species) + (1 | sp_day_year) + (1 | IDLocality),
data = d[Country == "Hungary"], REML = FALSE
)
est_hs <- est_out(hs, "Hungary: (1|weekday)+(1|genus)+(1|Species)+(1|sp_day_year)+(scale(Covid)|IDLocality)+(1|sp_loc)")
est_hs[, control_for_starting_distance := "yes"]
hx <- lmer(scale(log(FID)) ~
# scale(log(SD)) +
scale(log(FlockSize)) +
scale(log(BodyMass)) +
scale(sin(rad)) + scale(cos(rad)) +
# scale(Day)+
scale(Temp) +
scale(Covid) +
(1 | Year) + (1 | weekday) + (1 | genus) + (1 | Species) + (1 | sp_day_year) + (scale(Covid) | IDLocality) + (1 | sp_loc),
# (1 | Year) + (1 | weekday) + (1 | genus) + (1 | Species) + (1 | sp_day_year) + (1 | IDLocality),
data = d[Country == "Hungary"], REML = FALSE
)
est_hx <- est_out(hx, "Hungary: (1|weekday)+(1|genus)+(1|Species)+(1|sp_day_year)+(scale(Covid)|IDLocality)+(1|sp_loc)")
est_hx[, control_for_starting_distance := "no"]
# AU - singular fits only due to Year and random slope estimated as zero (removing those changes no results)
as <- lmer(scale(log(FID)) ~
scale(log(SD)) +
scale(log(FlockSize)) +
scale(log(BodyMass)) +
scale(sin(rad)) + scale(cos(rad)) +
# scale(Day)+
scale(Temp) +
scale(Covid) +
(1 | Year) +(1 | weekday) + (1 | genus) + (1 | Species) + (1 | sp_day_year) + (scale(Covid) | IDLocality) + (1 | sp_loc),
# (1 | Year) + (1 | weekday) + (1|genus) + (1 | Species) + (1 | sp_day_year) + (1 | IDLocality),
data = d[Country == "Australia"], REML = FALSE
)
est_as <- est_out(as, "Australia: (1|weekday)+(1|genus)+(1|Species)+(1|sp_day_year)+(scale(Covid)|IDLocality)+(1|sp_loc)")
est_as[, control_for_starting_distance := "yes"]
ax <- lmer(scale(log(FID)) ~
# scale(log(SD)) +
scale(log(FlockSize)) +
scale(log(BodyMass)) +
scale(sin(rad)) + scale(cos(rad)) +
# scale(Day)+
scale(Temp) +
scale(Covid) +
(1 | Year) + (1 | weekday) + (1 | genus) + (1 | Species) + (1 | sp_day_year) + (scale(Covid) | IDLocality) + (1 | sp_loc),
# (1 | Year) + (1 | weekday) + (1 | genus) + (1 | Species) + (1 | sp_day_year) + (1 | IDLocality),
data = d[Country == "Australia"], , REML = FALSE,
control = lmerControl(
optimizer ='optimx', optCtrl=list(method='nlminb'))
)
est_ax <- est_out(ax, "Australia: (1|weekday)+(1|genus)+(1|Species)+(1|sp_day_year)+(scale(Covid)|IDLocality)+(1|sp_loc)")
est_ax[, control_for_starting_distance := "no"]
# PL
ps <- lmer(scale(log(FID)) ~
scale(log(SD)) +
scale(log(FlockSize)) +
scale(log(BodyMass)) +
scale(sin(rad)) + scale(cos(rad)) +
# scale(Day)+
scale(Temp) +
scale(Covid) +
(1 | Year) + (1 | weekday) + (1 | genus) + (1 | Species) + (1 | sp_day_year),
# (1 | Year) + (1 | weekday) + (1|genus) + (1 | Species) + (1 | sp_day_year) + (1 | IDLocality),
data = d[Country == "Poland"], REML = FALSE
)
est_ps <- est_out(ps, "Poland: (1|weekday)+(1|genus)+(1|Species)+(1|sp_day_year)")
est_ps[, control_for_starting_distance := "yes"]
px <- lmer(scale(log(FID)) ~
# scale(log(SD)) +
scale(log(FlockSize)) +
scale(log(BodyMass)) +
scale(sin(rad)) + scale(cos(rad)) +
# scale(Day)+
scale(Temp) +
scale(Covid) +
(1 | Year) + (1 | weekday) + (1 | genus) + (1 | Species) + (1 | sp_day_year),
# (1 | Year) + (1 | weekday) + (1 | genus) + (1 | Species) + (1 | sp_day_year) + (1 | IDLocality),
data = d[Country == "Poland"], REML = FALSE
)
est_px <- est_out(px, "Poland: (1|weekday)+(1|genus)+(1|Species)+(1|sp_day_year)")
est_px[, control_for_starting_distance := "no"]
# combine
est_ms[, Country := 'All\n(mixed model)']
est_mx[, Country := "All\n(mixed model)"]
est_as[, Country := "Australia"]
est_ax[, Country := "Australia"]
est_cs[, Country := "Czechia"]
est_cx[, Country := "Czechia"]
est_hs[, Country := "Hungary"]
est_hx[, Country := "Hungary"]
est_ps[, Country := "Poland"]
est_px[, Country := "Poland"]
est_fs[, Country := "Finland"]
est_fx[, Country := "Finland"]
o = rbind(est_ms, est_mx,
est_as, est_ax,
est_cs, est_cx,
est_hs, est_hx,
est_ps, est_px,
est_fs, est_fx)
save(o, file = here::here('Data/dat_est_rev.Rdata'))
#AIC(hs, hx)
#AIC(cs, cx)
#0AIC(ps, px)
#AIC(as, ax)
#AIC(fs, fx)
# supplementary Sable S_fid_period
ms_out <- m_out(name = "Table S1 - full a", dep = "Escape distance", model = ms, nsim = 5000)
mx_out <- m_out(name = "Table S1 - full b", dep = "Escape distance", model = mx, nsim = 5000)
cs_out <- m_out(name = "Table S1 - CZ a", dep = "Escape distancey", model = cs, nsim = 5000)
cx_out <- m_out(name = "Table S1 - CZ b", dep = "Escape distance", model = cx, nsim = 5000)
fs_out <- m_out(name = "Table S1 - FI a", dep = "Escape distance", model = fs, nsim = 5000)
fx_out <- m_out(name = "Table S1 - FI b", dep = "Escape distance", model = fx, nsim = 5000)
hs_out <- m_out(name = "Table S1 - HU a", dep = "Escape distance", model = hs, nsim = 5000)
hx_out <- m_out(name = "Table S1 - HU b", dep = "Escape distance", model = hx, nsim = 5000)
as_out <- m_out(name = "Table S1 - AU a", dep = "Escape distance", model = as, nsim = 5000)
ax_out <- m_out(name = "Table S1 - AU b", dep = "Escape distance", model = ax, nsim = 5000)
ps_out <- m_out(name = "Table S1 - PL a", dep = "Escape distance", model = ps, nsim = 5000)
px_out <- m_out(name = "Table S1 - PL b", dep = "Escape distance", model = px, nsim = 5000)
out_FID_c <- rbind(fs_out, fx_out, ps_out, px_out, cs_out, cx_out, hs_out, hx_out, as_out, ax_out, fill = TRUE)
out_FID_c[is.na(out_FID_c)] <- ""
out_FID_c$R2_mar = out_FID_c$R2_con = NULL
out_FID_c[, effect := gsub("scale\\(Covid\\)", "Period", effect)]
out_FID_c[, effect := gsub("scale\\(Year\\)", "year", effect)]
out_FID_c[, effect := gsub("scale\\(log\\(SD\\)\\)", "starting distance (ln)", effect)]
out_FID_c[, effect := gsub("scale\\(Temp\\)", "temperaturre", effect)]
out_FID_c[, effect := gsub("scale\\(log\\(FlockSize\\)\\)", "flock size (ln)", effect)]
out_FID_c[, effect := gsub("scale\\(log\\(BodyMass\\)\\)", "body mass (ln)", effect)]
out_FID_c[, effect := gsub("scale\\(sin\\(rad\\)\\)", "time (sine of radians)", effect)]
out_FID_c[, effect := gsub("scale\\(cos\\(rad\\)\\)", "time (cosine of radians)", effect)]
out_FID_c[, effect := gsub("Species", "species", effect)]
out_FID_c[, effect := gsub("Year", "year", effect)]
out_FID_c[, effect := gsub("sp_day_year", "species within day & year", effect)]
out_FID_c[, effect := gsub("IDLocality", "site", effect)]
out_FID_c[, effect := gsub("sp_loc", "species within site", effect)]
out_FID_c[, effect := gsub("site Period", "Period (slope) | site", effect)]
fwrite(file = here::here("Outputs/Table_S1.csv"), out_FID_c)load(here::here("Data/dat_est_rev.Rdata"))
o[predictor %in% c("scale(Covid)"), predictor := "Period"]
oo <- o[predictor %in% c("Period")]
oo[, N:=as.numeric(sub('.*N = ', '', model))]
# add meta-analytical mean
oo_s = oo[control_for_starting_distance == 'yes']
met = summary(meta.summaries(d = oo_s$estimate, se = oo_s$sd, method = "fixed", weights = oo_s$N))$summci
oo_met = data.table(predictor = "Period", estimate = met[2], lwr = met[1], upr = met[3], sd = NA, model = NA, control_for_starting_distance = "yes", Country = "Combined\n(metanalytical)", N = NA)
oo_sx = oo[control_for_starting_distance == "no"]
metx = summary(meta.summaries(d = oo_sx$estimate, se = oo_sx$sd, method = "fixed", weights = oo_sx$N))$summci
oo_metx = data.table(predictor = "Period", estimate = metx[2], lwr = metx[1], upr = metx[3], sd = NA, model = NA, control_for_starting_distance = "no", Country = "Combined\n(metanalytical)", N = NA)
oo = rbind(oo, oo_met, oo_metx)
oo[, Country := factor(Country, levels = rev(c("Finland", "Poland", "Czechia", "Hungary", "Australia", "Combined\n(metanalytical)", "All\n(mixed model)")))]
# prepare for adding N
oo[control_for_starting_distance == "no" | is.na(N), N := ""]
oo[, n_pos := 1.1]
width_ <- .5 # spacing between error bars
#col_ <- c(brewer.pal(n = 12, name = "Paired"), "grey30", "grey80")
#Tol_bright <- c("#EE6677", "#228833", "#4477AA", "#CCBB44", "#66CCEE", "#AA3377", "#BBBBBB")
#Tol_muted <- c("#88CCEE", "#44AA99", "#117733", "#332288", "#DDCC77", "#999933", "#CC6677", "#882255", "#AA4499", "#DDDDDD")
#Tol_light <- c("#BBCC33", "#AAAA00", "#77AADD", "#EE8866", "#EEDD88", "#FFAABB", "#99DDFF", "#44BB99", "#DDDDDD")
# From Color Universal Design (CUD): https://jfly.uni-koeln.de/color/
#Okabe_Ito <- c("#E69F00", "#56B4E9", "#009E73", "#F0E442", "#0072B2", "#D55E00", "#CC79A7", "#000000")
#col_ = Okabe_Ito[7:1]
# JAMA and LocusZoom modified order
#col_ = c("#374E55FF", "#374E55FF", "#DF8F44FF", "#79AF97FF", "#00A1D5FF", "#B24745FF", "#80796BFF") #"#6A6599FF",
#col_ <- c("#357EBDFF", "#9632B8FF", "#46B8DAFF", "#5CB85CFF", "#EEA236FF", "#D43F3AFF", "#D43F3AFF")[7:1] # "#D43F3AFF", "#B8B8B8FF"
col_ = c("#357EBDFF", "#D43F3AFF", "#46B8DAFF", "#5CB85CFF", "#EEA236FF", "#9632B8FF", "#9632B8FF")[7:1] # "#D43F3AFF", "#B8B8B8FF"
#show_col(col_)
#p =
ggplot(oo, aes(x = estimate, y = Country, col = Country, shape = control_for_starting_distance)) +
geom_vline(xintercept = 0, color = "grey", linetype = "dotted") +
geom_errorbarh(aes(xmin = lwr, xmax = upr), height = 0, position = ggstance::position_dodgev(width_)) +
# geom_point(position = ggstance::position_dodgev(.6)) +
geom_point(position = position_dodge(width = width_), bg = "white", size = 1.1) +
# scale_color_viridis(discrete=TRUE, begin=0, end = 0.5) +
# scale_fill_viridis(discrete=TRUE, begin=0, end = 0.5) +
# geom_text( aes(x = n_pos,label = N), vjust = 0, size = 1.75, position = ggstance::position_dodgev(width_))+ # 3 positions for 3 bars
# annotate("text", x=log10(3), y=85, label= "Used", col = "grey30", size = 2.5)+
geom_text( aes(x = n_pos,label = N), vjust = 1, size = 1.75, position = ggstance::position_dodgev(width_))+
scale_shape_manual(name = "Controlled for\nstarting distance", guide = guide_legend(reverse = TRUE), values = c(21, 19)) +
#scale_color_jama(guide = "none")+ #, palette = 'light'
scale_color_manual(guide = "none", values = col_) + #guide_legend(reverse = TRUE)
scale_x_continuous(breaks = round(seq(-0.6, 1.2, by = 0.3), 1)) +
ylab("") +
xlab("Standardized effect size of Period\n[on flight initiation distance]") +
# coord_cartesian(xlim = c(-.15, .15)) +
# scale_x_continuous(breaks = round(seq(-.15, .15, by = 0.05),2)) +
theme_bw() +
theme(
legend.position = "right",
legend.title = element_text(size = 7),
legend.text = element_text(size = 6),
# legend.spacing.y = unit(0.1, 'cm'),
legend.key.height = unit(0.5, "line"),
legend.margin = margin(0, 0, 0, 0),
# legend.position=c(0.5,1.6),
plot.title = element_text(color = "grey", size = 7),
plot.margin = margin(b = 0.5, l = 0.5, t = 0.5, r = 0.5, unit = "pt"),
panel.grid = element_blank(),
panel.border = element_blank(),
panel.background = element_blank(),
axis.line = element_line(colour = ax_lines, size = 0.25),
axis.line.y = element_blank(),
axis.ticks.y = element_blank(),
axis.ticks.x = element_line(colour = ax_lines, size = 0.25),
# axis.text.x = element_text()
axis.ticks.length = unit(1, "pt"),
axis.text.x = element_text(, size = 6),
axis.text.y = element_text(colour = "black", size = 7),
axis.title = element_text(size = 7)
)
ggsave(here::here("Outputs/Fig_rev_width_CustomLocusZoom_v2.png"), width = 8, height = 6.5, unit = "cm", dpi = 600)
#ggsave("Outputs/Fig_rev_width_CustomJAMAv1.png", width = 8, height = 6, unit = "cm", dpi = 600)
#ggsave("Outputs/Fig_rev_width_Okabe_v2.png", width = 8, height = 6, unit = "cm", dpi = 600)
#ggsave("Outputs/Fig_rev_width_UChicago_v3.png", width = 8, height = 6, unit = "cm", dpi = 600)
#ggsave("Outputs/Fig_rev_width_d3_v2.png", width = 8, height = 6, unit = "cm", dpi = 600)
#ggsave("Outputs/Fig_rev_width_nejm_v2.png", width = 8, height = 6, unit = "cm", dpi = 600)
#ggsave("Outputs/Fig_rev_width_jama_v2.png", width = 8, height = 6, unit = "cm", dpi = 600)
#ggsave("Outputs/Fig_rev_width_jco_v2.png", width = 8, height = 6, unit = "cm", dpi = 600)
#ggsave("Outputs/Fig_rev_width_npg_v2.png", width = 8, height = 6, unit = "cm", dpi = 600)
# show used colors
#gg <- ggplot_build(p)
#col_ = unique(gg$data[[3]]["colour"])$colour
#show_col(col_)Figure 1 | Change in avian tolerance towards humans before vs during the COVID-19 shutdowns. The dots with horizontal lines represent estimated standardised effect size and their 95% confidence intervals, the numbers sample sizes. For the country-specific and "All“, the effect sizes and 95% confidence intervals come from the joint posterior distribution of 5000 simulated values generated by the sim function from the arm package (Gelman et al., 2022) using the mixed model outputs controlled for starting distance of the observer (filled circles) or not (empty circles; Tables S1 and S2a). The models were further controlled for flock size, body size, temperature (also a proxy for a day within the breeding season: rPearson = 0.48; Fig. S2), and time of a day, as well as for the non-independence of data points by fitting random intercepts of year, weekday, genus, species, species at a given day and year, country (in All - a global mixed model), site, and species within a site, while fitting Period as random slope within site (i.e. allowing for different Period effect at each site) and in All also within country. Fitting Period as random slope at other random intercepts produces similar results (see Fig. S1a). The multicollinearity was small as correlations between predictors were weak (Fig. S2). For the “Combined“, the estimate and 95% confidence interval represent the meta-analytical mean based on the country estimates and their standard deviation (from the country-specific models), and sample size per country. Note that effect sizes are small and estimates centre around zero.
Table S1 | Escape distance in relations to Period, given country
out_FID_c$error_structure = out_FID_c$response = NULL
out_FID_c[model!="", model:=c('Finland', 'Finland, without starting distance',
'Poland', 'Poland, without starting distance',
'Czechia', 'Czechia, without starting distance',
'Hungary', 'Hungary, without starting distance',
'Australia', 'Australia, without starting distance')]
setnames(out_FID_c, old = c("estimate_r", "lwr_r", "upr_r"), new = c("estimate", "lower", "upper"))
out_FID_c %>%
kbl() %>%
kable_paper("hover", full_width = F)| model | N | type | effect | estimate | lower | upper |
|---|---|---|---|---|---|---|
| Finland | 1019 | fixed | (Intercept) | -0.035 | -0.321 | 0.231 |
| fixed | starting distance (ln) | 0.282 | 0.228 | 0.336 | ||
| fixed | flock size (ln) | -0.049 | -0.101 | 0.001 | ||
| fixed | body mass (ln) | 0.091 | -0.066 | 0.24 | ||
| fixed | time (sine of radians) | 0.085 | 0.022 | 0.147 | ||
| fixed | time (cosine of radians) | 0.031 | -0.028 | 0.088 | ||
| fixed | temperaturre | -0.079 | -0.148 | -0.009 | ||
| fixed | Period | -0.15 | -0.337 | 0.031 | ||
| random | species within day & year (Intercept) | 10% | 9% | 11% | ||
| random | species within site (Intercept) | 6% | 5% | 7% | ||
| random | site (Intercept) | 1% | 0% | 3% | ||
| random | Period (slope) | site | 1% | 0% | 3% | ||
| random | species (Intercept) | 4% | 3% | 5% | ||
| random | genus (Intercept) | 10% | 7% | 14% | ||
| random | weekday (Intercept) | 0% | 0% | 1% | ||
| random | year (Intercept) | 5% | 1% | 11% | ||
| random | Residual | 61% | 50% | 72% | ||
| Finland, without starting distance | 1019 | fixed | (Intercept) | -0.125 | -0.468 | 0.221 |
| fixed | flock size (ln) | -0.018 | -0.07 | 0.034 | ||
| fixed | body mass (ln) | 0.2 | 0.03 | 0.375 | ||
| fixed | time (sine of radians) | 0.093 | 0.026 | 0.161 | ||
| fixed | time (cosine of radians) | 0.029 | -0.032 | 0.089 | ||
| fixed | temperaturre | -0.076 | -0.148 | -0.005 | ||
| fixed | Period | -0.132 | -0.366 | 0.101 | ||
| random | species within day & year (Intercept) | 9% | 7% | 10% | ||
| random | species within site (Intercept) | 7% | 6% | 7% | ||
| random | site (Intercept) | 1% | 1% | 2% | ||
| random | Period (slope) | site | 1% | 1% | 2% | ||
| random | species (Intercept) | 7% | 5% | 8% | ||
| random | genus (Intercept) | 10% | 7% | 13% | ||
| random | weekday (Intercept) | 1% | 0% | 1% | ||
| random | year (Intercept) | 9% | 2% | 16% | ||
| random | Residual | 56% | 44% | 68% | ||
| Poland | 762 | fixed | (Intercept) | 0.145 | -0.115 | 0.404 |
| fixed | starting distance (ln) | 0.479 | 0.419 | 0.538 | ||
| fixed | flock size (ln) | 0.012 | -0.042 | 0.066 | ||
| fixed | body mass (ln) | -0.07 | -0.172 | 0.031 | ||
| fixed | time (sine of radians) | -0.035 | -0.117 | 0.05 | ||
| fixed | time (cosine of radians) | -0.012 | -0.101 | 0.078 | ||
| fixed | temperaturre | -0.111 | -0.182 | -0.04 | ||
| fixed | Period | 0.119 | -0.112 | 0.353 | ||
| random | species within day & year (Intercept) | 8% | 8% | 8% | ||
| random | species (Intercept) | 15% | 12% | 17% | ||
| random | genus (Intercept) | 0% | 0% | 0% | ||
| random | weekday (Intercept) | 0% | 0% | 1% | ||
| random | year (Intercept) | 11% | 7% | 17% | ||
| random | Residual | 66% | 57% | 72% | ||
| Poland, without starting distance | 762 | fixed | (Intercept) | 0.246 | -0.075 | 0.558 |
| fixed | flock size (ln) | 0.054 | -0.009 | 0.117 | ||
| fixed | body mass (ln) | 0.046 | -0.097 | 0.19 | ||
| fixed | time (sine of radians) | -0.014 | -0.112 | 0.077 | ||
| fixed | time (cosine of radians) | -0.001 | -0.104 | 0.099 | ||
| fixed | temperaturre | -0.104 | -0.187 | -0.02 | ||
| fixed | Period | 0.245 | 0.001 | 0.498 | ||
| random | species within day & year (Intercept) | 11% | 10% | 11% | ||
| random | species (Intercept) | 26% | 22% | 29% | ||
| random | genus (Intercept) | 1% | 1% | 2% | ||
| random | weekday (Intercept) | 0% | 0% | 0% | ||
| random | year (Intercept) | 8% | 5% | 14% | ||
| random | Residual | 53% | 46% | 60% | ||
| Czechia | 2013 | fixed | (Intercept) | 0.111 | -0.257 | 0.481 |
| fixed | starting distance (ln) | 0.392 | 0.351 | 0.437 | ||
| fixed | flock size (ln) | 0.005 | -0.03 | 0.037 | ||
| fixed | body mass (ln) | 0.17 | 0.038 | 0.308 | ||
| fixed | time (sine of radians) | 0.05 | 0.002 | 0.099 | ||
| fixed | time (cosine of radians) | 0.054 | 0.009 | 0.097 | ||
| fixed | temperaturre | 0.043 | -0.008 | 0.093 | ||
| fixed | Period | 0.036 | -0.294 | 0.361 | ||
| random | species within day & year (Intercept) | 3% | 2% | 3% | ||
| random | species within site (Intercept) | 3% | 2% | 3% | ||
| random | species (Intercept) | 13% | 12% | 13% | ||
| random | site (Intercept) | 2% | 0% | 3% | ||
| random | Period (slope) | site | 2% | 0% | 3% | ||
| random | genus (Intercept) | 12% | 11% | 11% | ||
| random | weekday (Intercept) | 0% | 0% | 0% | ||
| random | year (Intercept) | 15% | 1% | 33% | ||
| random | Residual | 52% | 34% | 70% | ||
| Czechia, without starting distance | 2013 | fixed | (Intercept) | -0.043 | -0.969 | 0.867 |
| fixed | flock size (ln) | 0.005 | -0.031 | 0.04 | ||
| fixed | body mass (ln) | 0.231 | 0.078 | 0.382 | ||
| fixed | time (sine of radians) | 0.14 | 0.086 | 0.195 | ||
| fixed | time (cosine of radians) | 0.123 | 0.075 | 0.172 | ||
| fixed | temperaturre | 0.035 | -0.027 | 0.098 | ||
| fixed | Period | 0.223 | -0.642 | 1.1 | ||
| random | species within day & year (Intercept) | 2% | 1% | 5% | ||
| random | species within site (Intercept) | 1% | 1% | 2% | ||
| random | species (Intercept) | 6% | 4% | 10% | ||
| random | site (Intercept) | 1% | -1% | 1% | ||
| random | Period (slope) | site | 1% | -1% | 1% | ||
| random | genus (Intercept) | 8% | 5% | 11% | ||
| random | weekday (Intercept) | 1% | 0% | 1% | ||
| random | year (Intercept) | 54% | 17% | 73% | ||
| random | Residual | 26% | 13% | 57% | ||
| Hungary | 1456 | fixed | (Intercept) | 0.097 | -0.174 | 0.384 |
| fixed | starting distance (ln) | 0.485 | 0.438 | 0.531 | ||
| fixed | flock size (ln) | -0.012 | -0.051 | 0.031 | ||
| fixed | body mass (ln) | -0.065 | -0.184 | 0.059 | ||
| fixed | time (sine of radians) | -0.077 | -0.133 | -0.019 | ||
| fixed | time (cosine of radians) | -0.011 | -0.06 | 0.037 | ||
| fixed | temperaturre | -0.023 | -0.083 | 0.037 | ||
| fixed | Period | -0.091 | -0.271 | 0.097 | ||
| random | species within day & year (Intercept) | 5% | 4% | 5% | ||
| random | species within site (Intercept) | 0% | 0% | 0% | ||
| random | species (Intercept) | 3% | 2% | 3% | ||
| random | genus (Intercept) | 7% | 5% | 8% | ||
| random | site (Intercept) | 2% | 0% | 8% | ||
| random | Period (slope) | site | 2% | 0% | 8% | ||
| random | weekday (Intercept) | 1% | 0% | 1% | ||
| random | year (Intercept) | 5% | 2% | 8% | ||
| random | Residual | 76% | 59% | 84% | ||
| Hungary, without starting distance | 1456 | fixed | (Intercept) | 0.078 | -0.336 | 0.494 |
| fixed | flock size (ln) | -0.011 | -0.057 | 0.035 | ||
| fixed | body mass (ln) | 0.187 | 0.006 | 0.376 | ||
| fixed | time (sine of radians) | -0.111 | -0.18 | -0.044 | ||
| fixed | time (cosine of radians) | -0.007 | -0.06 | 0.05 | ||
| fixed | temperaturre | -0.025 | -0.095 | 0.043 | ||
| fixed | Period | -0.105 | -0.307 | 0.109 | ||
| random | species within day & year (Intercept) | 7% | 5% | 8% | ||
| random | species within site (Intercept) | 0% | 0% | 0% | ||
| random | species (Intercept) | 2% | 2% | 2% | ||
| random | genus (Intercept) | 18% | 13% | 18% | ||
| random | site (Intercept) | 1% | 0% | 12% | ||
| random | Period (slope) | site | 1% | 0% | 12% | ||
| random | weekday (Intercept) | 1% | 0% | 1% | ||
| random | year (Intercept) | 5% | 3% | 7% | ||
| random | Residual | 63% | 42% | 74% | ||
| Australia | 1119 | fixed | (Intercept) | -0.064 | -0.202 | 0.081 |
| fixed | starting distance (ln) | 0.502 | 0.452 | 0.553 | ||
| fixed | flock size (ln) | 0.02 | -0.026 | 0.068 | ||
| fixed | body mass (ln) | -0.066 | -0.163 | 0.029 | ||
| fixed | time (sine of radians) | -0.032 | -0.099 | 0.036 | ||
| fixed | time (cosine of radians) | -0.026 | -0.093 | 0.041 | ||
| fixed | temperaturre | 0.009 | -0.045 | 0.063 | ||
| fixed | Period | -0.009 | -0.08 | 0.062 | ||
| random | species within day & year (Intercept) | 7% | 6% | 7% | ||
| random | species within site (Intercept) | 13% | 11% | 13% | ||
| random | species (Intercept) | 14% | 11% | 14% | ||
| random | genus (Intercept) | 2% | 1% | 2% | ||
| random | site (Intercept) | 0% | -1% | 8% | ||
| random | Period (slope) | site | 0% | -1% | 8% | ||
| random | weekday (Intercept) | 0% | 0% | 0% | ||
| random | year (Intercept) | 0% | 0% | 0% | ||
| random | Residual | 64% | 51% | 68% | ||
| Australia, without starting distance | 1119 | fixed | (Intercept) | -0.18 | -0.432 | 0.058 |
| fixed | flock size (ln) | 0.051 | -0.002 | 0.104 | ||
| fixed | body mass (ln) | 0.055 | -0.073 | 0.184 | ||
| fixed | time (sine of radians) | -0.017 | -0.098 | 0.063 | ||
| fixed | time (cosine of radians) | -0.006 | -0.08 | 0.069 | ||
| fixed | temperaturre | 0.004 | -0.058 | 0.067 | ||
| fixed | Period | 0.156 | -0.007 | 0.331 | ||
| random | species within day & year (Intercept) | 8% | 5% | 10% | ||
| random | species within site (Intercept) | 12% | 8% | 15% | ||
| random | species (Intercept) | 17% | 14% | 17% | ||
| random | genus (Intercept) | 5% | 4% | 5% | ||
| random | site (Intercept) | 1% | -10% | 15% | ||
| random | Period (slope) | site | 1% | -10% | 15% | ||
| random | weekday (Intercept) | 0% | 0% | 0% | ||
| random | year (Intercept) | 1% | 0% | 2% | ||
| random | Residual | 55% | 36% | 73% |
For model descriptions see legend of Fig. 1, for descriptions of variables Methods of the paper.
Alternative models give similar results
# prepare estimates Period
# 01a all data, main text
m1a <- lmer(scale(log(FID)) ~
scale(log(SD)) +
scale(log(FlockSize)) +
scale(log(BodyMass)) +
scale(sin(rad)) + scale(cos(rad)) +
# scale(Day)+
scale(Temp) +
scale(Covid) +
(1 | Year) + (1 | weekday) + (1 | genus) + (1 | Species) + (1 | sp_day_year) + (scale(Covid) | Country) + (scale(Covid) | IDLocality) + (1 | sp_loc),
data = d, REML = FALSE,
control = lmerControl(
optimizer = "optimx", optCtrl = list(method = "nlminb")
)
)
est_m1a = est_out(m1a, "01a) (1|Year) + (1|weekday) + (1|genus) + (1|Species) + (1|sp_day_year) + (scale(Covid)|Country) + (scale(Covid)|IDLocality) + (1|sp_loc)")
# 01b all data, all random slopes - singularity
m1b=lmer(scale(log(FID))~
scale(log(SD))+
scale(log(FlockSize))+
scale(log(BodyMass))+
scale(sin(rad)) + scale(cos(rad)) +
#scale(Day)+
scale(Temp)+
scale(Covid)+
(1|Year) + (1 | weekday) + (scale(Covid)|genus)+(scale(Covid)|Species)+(1|sp_day_year) + (scale(Covid)|Country) +(scale(Covid)|IDLocality) + (scale(Covid)|sp_loc),
data = d, REML = FALSE
)
# # (Covid|IDLocality) +
est_m1b = est_out(m1b, "01b) (1|Year) + (1|weekday) + (scale(Covid)|genus) + (scale(Covid)|Species) + (1|sp_day_year) + (scale(Covid)|Country) + (scale(Covid)|IDLocality) + (scale(Covid)|sp_loc)")
# 01c all data, all random slopes, but some without cor to avoid singularity
m1c=lmer(scale(log(FID))~
scale(log(SD))+
scale(log(FlockSize))+
scale(log(BodyMass))+
scale(sin(rad)) + scale(cos(rad)) +
#scale(Day)+
scale(Temp)+
scale(Covid)+
#(1|Year) +(0+Covid|genus)+(0+Covid|Species)+(1|sp_day_year) + (Covid|Country) + (0+Covid|IDLocality) +(Covid|sp_loc)
(1|Year)+ (1|weekday) + (0+scale(Covid)|genus)+(0+scale(Covid)|Species)+(1|sp_day_year) + (scale(Covid)|Country) +(scale(Covid)|IDLocality) + (scale(Covid)|sp_loc),
data = d, REML = FALSE) # (0+scale(Covid)|genus) is zero
est_m1c = est_out(m1c, "01c) (1|Year) + (1|weekday) + (0+scale(Covid)|genus) + (0+scale(Covid)|Species) + (1|sp_day_year) + (scale(Covid)|Country) + (scale(Covid)|IDLocality) + (scale(Covid)|sp_loc)")
# 01d all data, random slopes that allow for non-singular fit
m1d=lmer(scale(log(FID))~
scale(log(SD))+
scale(log(FlockSize))+
scale(log(BodyMass))+
scale(sin(rad)) + scale(cos(rad)) +
#scale(Day)+
scale(Temp)+
scale(Covid)+
(1|Year) + (1|weekday) +(1|genus)+(1|Species)+(1|sp_day_year) + (scale(Covid)|Country) + (scale(Covid)|IDLocality) + (scale(Covid)|sp_loc),
data = d, REML = FALSE)# (Covid|IDLocality) +
#d[,res := resid(mf)]
est_m1d = est_out(m1d, "01d) (1|Year) + (1|weekday)+ (1|genus) + (1|Species) + (1|sp_day_year) + (scale(Covid)|Country) + (scale(Covid)|IDLocality) + (scale(Covid)|sp_loc)")
# 01e all data, random slopes that allow for non-singular fit with simple structure
m1e=lmer(scale(log(FID))~
scale(log(SD))+
scale(log(FlockSize))+
scale(log(BodyMass))+
scale(sin(rad)) + scale(cos(rad)) +
#scale(Day)+
scale(Temp)+
scale(Covid)+
(1|Year) + (1|weekday)+(1|genus)+(1|Species)+(1|sp_day_year) + (scale(Covid)|Country) + (1|IDLocality) + (1|sp_loc),
data = d, REML = FALSE)# (Covid|IDLocality) +
est_m1e = est_out(m1e, "01e) (1|Year) + (1|weekday) + (1|genus) + (1|Species) + (1|sp_day_year) + (scale(Covid)|Country) + (1|IDLocality) + (1|sp_loc)")
# 02a) before & during > 4/species - main text
dx = dd[N_during>4 & N_before >4]
#dxx = d[Species %in% dx$Species]
#dx2 = dd[N_during>9 & N_before >9]
m2a <- lmer(scale(log(FID)) ~
scale(log(SD)) +
scale(log(FlockSize)) +
scale(log(BodyMass)) +
scale(sin(rad)) + scale(cos(rad)) +
# scale(Day)+
scale(Temp) +
scale(Covid) +
(1 | Year) + (1| weekday) + (1 | genus) + (1 | Species) + (1 | sp_day_year) + (scale(Covid) | Country) + (scale(Covid) | IDLocality) + (1 | sp_loc),
data = d[Species %in% dx$Species], REML = FALSE,
control = lmerControl(
optimizer = "optimx", optCtrl = list(method = "nlminb")
)
)
est_m2a = est_out(m2a, "02a) (1|Year) + (1|weekday) + (1|genus) + (1|Species) + (1|sp_day_year) + (scale(Covid)|Country) + (scale(Covid)|IDLocality) + (1|sp_loc); >4/species/period")
# 02b) before & during > 4/species - singularity
dx = dd[N_during > 4 & N_before > 4]
m2b=lmer(scale(log(FID))~
scale(log(SD))+
scale(log(FlockSize))+
scale(log(BodyMass))+
scale(sin(rad)) + scale(cos(rad)) +
#scale(Day)+
scale(Temp)+
scale(Covid)+
(1|Year) + (1| weekday)+ (scale(Covid)|genus)+(scale(Covid)|Species)+(1|sp_day_year) + (scale(Covid)|Country) +(scale(Covid)|IDLocality) + (scale(Covid)|sp_loc),
#(1|Year) + (1|genus)+(1|Species)+(1|sp_day_year) + (Covid|Country) + (Covid|sp_loc),
data = d[Species %in% dx$Species],
REML = FALSE) # (Covid|IDLocality) +
est_m2b = est_out(m2b, "02b) (1|Year) + (1|weekday)+ (scale(Covid)|genus)+(scale(Covid)|Species) + (1|sp_day_year) + (scale(Covid)|Country) + (scale(Covid)|IDLocality) + (scale(Covid)|sp_loc); >4/species/period")
# 02c) before & during > 4/species, random slopes that allow for non-singular fit and simple structure
dx = dd[N_during>4 & N_before >4]
#dxx = d[Species %in% dx$Species]
#dx2 = dd[N_during>9 & N_before >9]
m2c=lmer(scale(log(FID))~
scale(log(SD))+
scale(log(FlockSize))+
scale(log(BodyMass))+
scale(sin(rad)) + scale(cos(rad)) +
#scale(Day)+
scale(Temp)+
scale(Covid)+
(1|Year)+ (1| weekday) +(1|genus)+(1|Species)+(1|sp_day_year) + (scale(Covid)|Country) + (1|IDLocality) + (1|sp_loc),
#(1|Year) + (1|genus)+(1|Species)+(1|sp_day_year) + (Covid|Country) + (Covid|sp_loc),
data = d[Species %in% dx$Species],
REML = FALSE) # (Covid|IDLocality) +
est_m2c = est_out(m2c, "02c) (1|Year) + (1|weekday)+ (1|genus) + (1|Species)+(1|sp_day_year) + (scale(Covid)|Country) + (1|IDLocality) + (1|sp_loc); >4/species/period")
# 03a) before & during > 9/species - main text
dx = dd[N_during > 9 & N_before > 9]
m3a <- lmer(scale(log(FID)) ~
scale(log(SD)) +
scale(log(FlockSize)) +
scale(log(BodyMass)) +
scale(sin(rad)) + scale(cos(rad)) +
# scale(Day)+
scale(Temp) +
scale(Covid) +
(1 | Year) + (1 | weekday) + (1 | genus) + (1 | Species) + (1 | sp_day_year) + (scale(Covid) | Country) + (scale(Covid) | IDLocality) + (1 | sp_loc),
data = d[Species %in% dx$Species], REML = FALSE,
)
est_m3a = est_out(m3a, "03a) (1|Year) + (1|weekday) + (1|genus) + (1|Species) + (1|sp_day_year) + (scale(Covid)|Country) + (scale(Covid)|IDLocality) + (1|sp_loc); >9/species/period")
# 03b) before & during > 9/species - singularity
dx = dd[N_during>9 & N_before >9]
#dxx = d[Species %in% dx$Species]
#dx2 = dd[N_during>9 & N_before >9]
m3b=lmer(scale(log(FID))~
scale(log(SD))+
scale(log(FlockSize))+
scale(log(BodyMass))+
scale(sin(rad)) + scale(cos(rad)) +
#scale(Day)+
scale(Temp)+
scale(Covid)+
(1|Year) + (1| weekday) + (scale(Covid)|genus)+(scale(Covid)|Species)+(1|sp_day_year) + (scale(Covid)|Country) +(scale(Covid)|IDLocality) + (scale(Covid)|sp_loc),
#(1|Year) + (1|genus)+(1|Species)+(1|sp_day_year) + (Covid|Country) + (Covid|sp_loc),
data = d[Species %in% dx$Species],
REML = FALSE) # (Covid|IDLocality) +
est_m3b = est_out(m3b, "03b) (1|Year) + (1|weekday) + (scale(Covid)|genus)+(scale(Covid)|Species) + (1|sp_day_year) + (scale(Covid)|Country) + (scale(Covid)|IDLocality) + (scale(Covid)|sp_loc); >9/species/period")
# 03c) before & during > 9/species, random slopes that allow for non-singular fit and simple structure
dx = dd[N_during>9 & N_before >9]
#dxx = d[Species %in% dx$Species]
#dx2 = dd[N_during>9 & N_before >9]
m3c=lmer(scale(log(FID))~
scale(log(SD))+
scale(log(FlockSize))+
scale(log(BodyMass))+
scale(sin(rad)) + scale(cos(rad)) +
#scale(Day)+
scale(Temp)+
scale(Covid)+
(1|Year) + (1| weekday) +(1|genus)+(1|Species)+(1|sp_day_year) + (scale(Covid)|Country) + (1|IDLocality) + (1|sp_loc),
#(1|Year) + (1|genus)+(1|Species)+(1|sp_day_year) + (Covid|Country) + (Covid|sp_loc),
data = d[Species %in% dx$Species],
REML = FALSE) # (Covid|IDLocality) +
est_m3c = est_out(m3c, '03c) (1|Year) + (1|genus) + (1|Species)+(1|sp_day_year) + (scale(Covid)|Country) + (1|IDLocality) + (1|sp_loc); >9/species/period')
# prepare estimates Stringency
m01a <- lmer(scale(log(FID)) ~
scale(Year) +
scale(log(SD)) +
scale(log(FlockSize)) +
scale(log(BodyMass)) +
scale(sin(rad)) + scale(cos(rad)) +
# scale(Day)+
scale(Temp) +
scale(StringencyIndex) +
(1|weekday) + (1 | genus) + (1 | Species) + (1 | sp_day_year) + (scale(StringencyIndex) | Country) + (1 | IDLocality) + (1 | sp_loc),
data = s, REML = FALSE)
est_m01a = est_out(m01a, "01a) (1|weekday) + (1|genus) + (1|Species) + (1|sp_day_year) + (scale(StringencyIndex)|Country) + (1|IDLocality) +(1|sp_loc)")
m01b=lmer(scale(log(FID))~
scale(Year)+
scale(log(SD))+
scale(log(FlockSize))+
scale(log(BodyMass))+
scale(sin(rad)) + scale(cos(rad)) +
#scale(Day)+
scale(Temp)+
scale(StringencyIndex)+
(1|weekday) + (scale(StringencyIndex)|genus)+(1|Species)+(1|sp_day_year) + (scale(StringencyIndex)|Country) + (1|IDLocality) +(1|sp_loc),
data = s, REML = FALSE
)
# (1|Year) explains nothing - could stay
est_m01b = est_out(m01b, "01b) (1|weekday) + (scale(StringencyIndex)|genus) + (1|Species) + (1|sp_day_year) + (scale(StringencyIndex)|Country) + (1|IDLocality) +(1|sp_loc)")
m01c=lmer(scale(log(FID))~
scale(Year)+
scale(log(SD))+
scale(log(FlockSize))+
scale(log(BodyMass))+
scale(sin(rad)) + scale(cos(rad)) +
#scale(Day)+
scale(Temp)+
scale(StringencyIndex)+
(scale(StringencyIndex)|Country) + (1|IDLocality),
data = s, REML = FALSE)
# (1|Year), (1|genus), (1|sp_day_year), (1|sp_loc), explain nothing - could stay
est_m01c = est_out(m01c, "01c) (scale(StringencyIndex)|Country) + (1|IDLocality)")
m02a <- lmer(scale(log(FID)) ~
scale(Year) +
scale(log(SD)) +
scale(log(FlockSize)) +
scale(log(BodyMass)) +
scale(sin(rad)) + scale(cos(rad)) +
# scale(Day)+
scale(Temp) +
scale(StringencyIndex) +
(1 | weekday) + (1 | genus) + (1 | Species) + (1 | sp_day_year) + (scale(StringencyIndex) | Country) + (1 | IDLocality) + (1 | sp_loc),
data = s[Nsp>4], REML = FALSE
)
est_m02a = est_out(m02a, "02a) (1|weekday) + (1|genus) + (1|Species) + (1|sp_day_year) + (scale(StringencyIndex)|Country) + (1|IDLocality) +(1|sp_loc); >4/species")
m02b=lmer(scale(log(FID))~
scale(Year)+
scale(log(SD))+
scale(log(FlockSize))+
scale(log(BodyMass))+
scale(sin(rad)) + scale(cos(rad)) +
#scale(Day)+
scale(Temp)+
scale(StringencyIndex)+
(1|weekday) + (scale(StringencyIndex)|genus)+(1|Species)+(1|sp_day_year) + (scale(StringencyIndex)|Country) + (1|IDLocality) +(1|sp_loc),
data = s[Nsp>4], REML = FALSE,
control = lmerControl(
optimizer ='optimx', optCtrl=list(method='nlminb'))
)
# (1|Year) explains nothing - could stay
est_m02b = est_out(m02b, "02b) (1|weekday) + (scale(StringencyIndex)|genus) + (1|Species) + (1|sp_day_year) + (scale(StringencyIndex)|Country) + (1|IDLocality) +(1|sp_loc); >4/species")
m02c=lmer(scale(log(FID))~
scale(Year)+
scale(log(SD))+
scale(log(FlockSize))+
scale(log(BodyMass))+
scale(sin(rad)) + scale(cos(rad)) +
#scale(Day)+
scale(Temp)+
scale(StringencyIndex)+
(scale(StringencyIndex)|Country) + (1|IDLocality),
data = s[Nsp>4], REML = FALSE)
# (1|Year), (1|genus), (1|sp_day_year), (1|sp_loc), explain nothing - could stay
est_m02c = est_out(m02c, "02c) (scale(StringencyIndex)|Country) + (1|IDLocality); >4/species")
m03a <- lmer(scale(log(FID)) ~
scale(Year) +
scale(log(SD)) +
scale(log(FlockSize)) +
scale(log(BodyMass)) +
scale(sin(rad)) + scale(cos(rad)) +
# scale(Day)+
scale(Temp) +
scale(StringencyIndex) +
(1 | weekday) + (1 | genus) + (1 | Species) + (1 | sp_day_year) + (scale(StringencyIndex) | Country) + (1 | IDLocality) + (1 | sp_loc),
data = s[Nsp > 9], REML = FALSE
)
est_m03a = est_out(m03a, "03a) (1|weekday) + (1|genus) + (1|Species) + (1|sp_day_year) + (scale(StringencyIndex)|Country) + (1|IDLocality) +(1|sp_loc); >9/species")
m03b = lmer(scale(log(FID)) ~
scale(Year) +
scale(log(SD)) +
scale(log(FlockSize)) +
scale(log(BodyMass)) +
scale(sin(rad)) + scale(cos(rad)) +
# scale(Day)+
scale(Temp) +
scale(StringencyIndex) +
(1 | weekday) + (scale(StringencyIndex) | genus) + (1 | Species) + (1 | sp_day_year) + (scale(StringencyIndex) | Country) + (1 | IDLocality) + (1 | sp_loc),
data = s[Nsp > 9], REML = FALSE
)
est_m03b = est_out(m03b, "03b) (1|weekday) + (scale(StringencyIndex)|genus) + (1|Species) + (1|sp_day_year) + (scale(StringencyIndex)|Country) + (1|IDLocality) +(1|sp_loc); >9/species")
m03c = lmer(scale(log(FID)) ~
scale(Year) +
scale(log(SD)) +
scale(log(FlockSize)) +
scale(log(BodyMass)) +
scale(sin(rad)) + scale(cos(rad)) +
# scale(Day)+
scale(Temp) +
scale(StringencyIndex) +
(scale(StringencyIndex) | Country) + (1 | IDLocality),
data = s[Nsp > 9], REML = FALSE
)
# (1|Year), (1|genus), (1|sp_day_year), (1|sp_loc), explain nothing - could stay
est_m03c = est_out(m03c, "03c) (scale(StringencyIndex)|Country) + 1|IDLocality); >9/species")
# prepare estimates google mobility
mg01a=lmer(scale(log(FID))~
scale(Year)+
scale(log(SD))+
scale(log(FlockSize))+
scale(log(BodyMass))+
scale(sin(rad)) + scale(cos(rad)) +
#scale(Day)+
scale(Temp)+
scale(parks_percent_change_from_baseline)+
(1|weekday) + (1|genus) +(1|Species) + (1|sp_day_year) + (scale(parks_percent_change_from_baseline)|Country) + (1|IDLocality)+(1|sp_loc),
data = ss, REML = FALSE
)
# (1|Year), (1|genus), (1|sp_day_year), (1|sp_loc), explain nothing - could stay
est_mg01a = est_out(mg01a, "01a) (1|weekday) + (1|genus) +(1|Species) + (1|sp_day_year) + (scale(parks_percent_change_from_baseline)|Country) + (1|IDLocality)+(1|sp_loc)")
mg01b <- lmer(scale(log(FID)) ~
scale(Year) +
scale(log(SD)) +
scale(log(FlockSize)) +
scale(log(BodyMass)) +
scale(sin(rad)) + scale(cos(rad)) +
# scale(Day)+
scale(Temp) +
scale(parks_percent_change_from_baseline) +
(1|weekday) + (scale(parks_percent_change_from_baseline)| genus) + (1 | Species) + (1 | sp_day_year) +
(scale(parks_percent_change_from_baseline)| Country) + (1 | IDLocality) + (1 | sp_loc),
data = ss, REML = FALSE
)
est_mg01b = est_out(mg01b, "01b) (1|weeekday) + (scale(parks_percent_change_from_baseline)|genus)+(1|Species)+(1|sp_day_year) + (scale(parks_percent_change_from_baseline)|Country) + (1|IDLocality) +(1|sp_loc)")
mg01c <- lmer(scale(log(FID)) ~
scale(Year) +
scale(log(SD)) +
scale(log(FlockSize)) +
scale(log(BodyMass)) +
scale(sin(rad)) + scale(cos(rad)) +
# scale(Day)+
scale(Temp) +
scale(parks_percent_change_from_baseline) +
(scale(parks_percent_change_from_baseline) | Country) + (1 | IDLocality),
data = ss, REML = FALSE
)
est_mg01c = est_out(mg01c, "01c) (scale(parks_percent_change_from_baseline)|Country) + (1|IDLocality)")
mg02a = lmer(scale(log(FID)) ~
scale(Year) +
scale(log(SD)) +
scale(log(FlockSize)) +
scale(log(BodyMass)) +
scale(sin(rad)) + scale(cos(rad)) +
# scale(Day)+
scale(Temp) +
scale(parks_percent_change_from_baseline) +
(1 | weekday) + (1 | genus) + (1 | Species) + (1 | sp_day_year) + (scale(parks_percent_change_from_baseline) | Country) + (1 | IDLocality) + (1 | sp_loc),
data = ss[Nsp>4], REML = FALSE
)
# (1|Year), (1|genus), (1|sp_day_year), (1|sp_loc), explain nothing - could stay
est_mg02a = est_out(mg02a, "02a) (1|weekday) + (1|genus) +(1|Species) + (1|sp_day_year) + (scale(parks_percent_change_from_baseline)|Country) + (1|IDLocality)+(1|sp_loc); >4/specie")
mg02b <- lmer(scale(log(FID)) ~
scale(Year) +
scale(log(SD)) +
scale(log(FlockSize)) +
scale(log(BodyMass)) +
scale(sin(rad)) + scale(cos(rad)) +
# scale(Day)+
scale(Temp) +
scale(parks_percent_change_from_baseline) +
(1 | weekday) + (scale(parks_percent_change_from_baseline) | genus) + (1 | Species) + (1 | sp_day_year) +
(scale(parks_percent_change_from_baseline) | Country) + (1 | IDLocality) + (1 | sp_loc),
data = ss[Nsp>4], REML = FALSE
)
est_mg02b = est_out(mg02b, "02b) (1|weeekday) + (scale(parks_percent_change_from_baseline)|genus)+(1|Species)+(1|sp_day_year) + (scale(parks_percent_change_from_baseline)|Country) + (1|IDLocality) +(1|sp_loc); >4/specie")
mg02c <- lmer(scale(log(FID)) ~
scale(Year) +
scale(log(SD)) +
scale(log(FlockSize)) +
scale(log(BodyMass)) +
scale(sin(rad)) + scale(cos(rad)) +
# scale(Day)+
scale(Temp) +
scale(parks_percent_change_from_baseline) +
(scale(parks_percent_change_from_baseline) | Country) + (1 | IDLocality),
data = ss[Nsp>4], REML = FALSE
)
est_mg02c = est_out(mg02c, "02c) (scale(parks_percent_change_from_baseline)|Country) + (1|IDLocality); >4/specie")
mg03a = lmer(scale(log(FID)) ~
scale(Year) +
scale(log(SD)) +
scale(log(FlockSize)) +
scale(log(BodyMass)) +
scale(sin(rad)) + scale(cos(rad)) +
# scale(Day)+
scale(Temp) +
scale(parks_percent_change_from_baseline) +
(1 | weekday) + (1 | genus) + (1 | Species) + (1 | sp_day_year) + (scale(parks_percent_change_from_baseline) | Country) + (1 | IDLocality) + (1 | sp_loc),
data = ss[Nsp > 9], REML = FALSE
)
# (1|Year), (1|genus), (1|sp_day_year), (1|sp_loc), explain nothing - could stay
est_mg03a = est_out(mg03a, "03a) (1|weekday) + (1|genus) +(1|Species) + (1|sp_day_year) + (scale(parks_percent_change_from_baseline)|Country) + (1|IDLocality)+(1|sp_loc); >9/specie")
mg03b <- lmer(scale(log(FID)) ~
scale(Year) +
scale(log(SD)) +
scale(log(FlockSize)) +
scale(log(BodyMass)) +
scale(sin(rad)) + scale(cos(rad)) +
# scale(Day)+
scale(Temp) +
scale(parks_percent_change_from_baseline) +
(1 | weekday) + (scale(parks_percent_change_from_baseline) | genus) + (1 | Species) + (1 | sp_day_year) +
(scale(parks_percent_change_from_baseline) | Country) + (1 | IDLocality) + (1 | sp_loc),
data = ss[Nsp > 9], REML = FALSE
)
est_mg03b = est_out(mg03b, "03b) (1|weeekday) + (scale(parks_percent_change_from_baseline)|genus)+(1|Species)+(1|sp_day_year) + (scale(parks_percent_change_from_baseline)|Country) + (1|IDLocality) +(1|sp_loc); >9/specie")
mg03c <- lmer(scale(log(FID)) ~
scale(Year) +
scale(log(SD)) +
scale(log(FlockSize)) +
scale(log(BodyMass)) +
scale(sin(rad)) + scale(cos(rad)) +
# scale(Day)+
scale(Temp) +
scale(parks_percent_change_from_baseline) +
(scale(parks_percent_change_from_baseline) | Country) + (1 | IDLocality),
data = ss[Nsp > 9], REML = FALSE
)
est_mg03c = est_out(mg03c, "03c) (scale(parks_percent_change_from_baseline)|Country) + (1|IDLocality); >9/specie")
# export
save(file = here::here("Data/Fig_S2_estimates.Rdata"),
est_m1a, est_m1b, est_m1c, est_m1d, est_m1e, est_m2a, est_m2b, est_m2c, est_m3a, est_m3b, est_m3c,
est_m01a, est_m01b, est_m01c, est_m02a, est_m02b, est_m02c, est_m03a, est_m03b, est_m03c,
est_mg01a, est_mg01b, est_mg01c, est_mg02a, est_mg02b, est_mg02c, est_mg03a, est_mg03b, est_mg03c)
# prepare plot for Period
xs = rbind(est_m1a, est_m1b, est_m1c, est_m1d, est_m1e, est_m2a, est_m2b, est_m2c, est_m3a, est_m3b, est_m3c)
#gsub("scale\\(Covid\\)", "Period", "(scale(Covid)|Country)")
xs[, model := gsub("scale\\(Covid\\)", "Period", xs$model)]
xs[, model := gsub("Year", "year", xs$model)]
xs[, model := gsub("Species", "species", xs$model)]
xs[, model := gsub("sp_day_year", "species within day & year", xs$model)]
xs[, model := gsub("IDLocality", "site", xs$model)]
xs[, model := gsub("sp_loc", "species within site", xs$model)]
gs2_ =
ggplot(xs[predictor == "scale(Covid)"], aes(y = model, x = estimate, col = model)) +
geom_vline(xintercept = 0, col = "grey30", lty = 3) +
geom_errorbar(aes(xmin = lwr, xmax = upr, col = model), width = 0, position = position_dodge(width = 0.01)) +
# ggtitle ("Sim based")+
geom_point(position = position_dodge(width = 0.01)) +
# scale_colour_brewer(type = 'qual', palette = 'Paired',guide = guide_legend(reverse = TRUE))+
# scale_fill_brewer(type = 'qual', palette = 'Paired',guide = guide_legend(reverse = TRUE))+
scale_color_viridis(discrete = TRUE, guide = guide_legend(reverse = TRUE)) +
scale_fill_viridis(discrete = TRUE, guide = guide_legend(reverse = TRUE)) +
scale_y_discrete(limits = rev) +
coord_fixed(ratio = 0.05, xlim = c(-0.23, 0.15)) +
# scale_shape(guide = guide_legend(reverse = TRUE)) +
# scale_x_continuous(limits = c(-2, 2), expand = c(0, 0), breaks = seq(-2,2, by = 1), labels = seq(-2,2, by = 1)) +
labs(y = NULL, x = "Period (before vs during shutdowns)", tag = "a)") + # title = "a)",Effect of Period (before/during shutdown)
# ylim(c(0,100))+
# coord_flip()+
theme_bw() +
theme(
legend.position = "none",
#plot.title.position = "plot",
plot.subtitle = element_text(size = 7),
plot.title = element_text(size = 7),
plot.tag = element_text(size = 7),
legend.title = element_text(size = 7),
legend.text = element_text(size = 6),
## legend.spacing.y = unit(0.1, 'cm'),
legend.key.height = unit(0.5, "line"),
# plot.margin = margin(b = 0.5, l = 0.5, t = 0.5, r =0.5, unit = "pt"),
panel.grid = element_blank(),
panel.border = element_blank(),
panel.background = element_blank(),
axis.line = element_line(colour = ax_lines, size = 0.25),
axis.line.y = element_blank(),
axis.ticks.y = element_blank(),
axis.ticks.x = element_line(colour = ax_lines, size = 0.25),
axis.ticks.length = unit(1, "pt"),
axis.text.x = element_text(colour = "black", size = 6),
axis.text.y = element_text(colour = "black", size = 7),
axis.title = element_text(size = 7)
)
#gs2
# ggsave(here::here('Outputs/Figure_Sy.png'),g, width = 30, height =5, units = 'cm')
# prepare plot for Stringency
xs0 = rbind(est_m01a, est_m01b, est_m01c, est_m02a, est_m02b, est_m02c, est_m03a, est_m03b, est_m03c)
xs0[, model := gsub("scale\\(StringencyIndex\\)", "Stringency Index", model)]
xs0[, model := gsub("Year", "year", model)]
xs0[, model := gsub("Species", "species", model)]
xs0[, model := gsub("sp_day_year", "species within day & year", model)]
xs0[, model := gsub("IDLocality", "site", model)]
xs0[, model := gsub("sp_loc", "species within site", model)]
g0_ =
ggplot(xs0[predictor == "scale(StringencyIndex)"], aes(y = model, x = estimate, col = model)) +
geom_vline(xintercept = 0, col = "grey30", lty = 3) +
geom_errorbar(aes(xmin = lwr, xmax = upr, col = model), width = 0, position = position_dodge(width = 0.01)) +
# ggtitle ("Sim based")+
geom_point(position = position_dodge(width = 0.01)) +
# scale_colour_brewer(type = 'qual', palette = 'Paired',guide = guide_legend(reverse = TRUE))+
# scale_fill_brewer(type = 'qual', palette = 'Paired',guide = guide_legend(reverse = TRUE))+
scale_color_viridis(discrete = TRUE, guide = guide_legend(reverse = TRUE)) +
scale_fill_viridis(discrete = TRUE, guide = guide_legend(reverse = TRUE)) +
scale_y_discrete(limits = rev) +
coord_fixed(ratio = 0.05, xlim = c(-0.23, 0.15)) +
# scale_shape(guide = guide_legend(reverse = TRUE)) +
# scale_x_continuous(limits = c(-2, 2), expand = c(0, 0), breaks = seq(-2,2, by = 1), labels = seq(-2,2, by = 1)) +
labs(y = NULL, x = "Stringency index", tag = "b)") + # title = "b) Effect of ") +
# ylim(c(0,100))+
# coord_flip()+
theme_bw() +
theme(
legend.position = "none",
plot.subtitle = element_text(size = 7),
plot.title = element_text(size = 7),
plot.tag = element_text(size = 7),
legend.title = element_text(size = 7),
legend.text = element_text(size = 6),
## legend.spacing.y = unit(0.1, 'cm'),
legend.key.height = unit(0.5, "line"),
# plot.margin = margin(b = 0.5, l = 0.5, t = 0.5, r =0.5, unit = "pt"),
panel.grid = element_blank(),
panel.border = element_blank(),
panel.background = element_blank(),
axis.line = element_line(colour = ax_lines, size = 0.25),
axis.line.y = element_blank(),
axis.ticks.y = element_blank(),
axis.ticks.x = element_line(colour = ax_lines, size = 0.25),
axis.ticks.length = unit(1, "pt"),
axis.text.x = element_text(colour = "black", size = 6),
axis.text.y = element_text(colour = "black", size = 7),
axis.title = element_text(size = 7)
)
#g0
# ggsave(here::here('Outputs/Figure_Sz.png'),g0, width = 30, height =5, units = 'cm')
# prepare plot for Google
xg0 = rbind(est_mg01a, est_mg01b, est_mg01c, est_mg02a, est_mg02b, est_mg02c, est_mg03a, est_mg03b, est_mg03c)
xg0[, model := gsub("scale\\(parks_percent_change_from_baseline\\)", "Google Mobility", model)]
xg0[, model := gsub("Year", "year", model)]
xg0[, model := gsub("Species", "species", model)]
xg0[, model := gsub("sp_day_year", "species within day & year", model)]
xg0[, model := gsub("IDLocality", "site", model)]
xg0[, model := gsub("sp_loc", "species within site", model)]
gg0_ =
ggplot(xg0[predictor == "scale(parks_percent_change_from_baseline)"], aes(y = model, x = estimate, col = model)) +
geom_vline(xintercept = 0, col = "grey30", lty = 3) +
geom_errorbar(aes(xmin = lwr, xmax = upr, col = model), width = 0, position = position_dodge(width = 0.01)) +
# ggtitle ("Sim based")+
geom_point(position = position_dodge(width = 0.01)) +
# scale_colour_brewer(type = 'qual', palette = 'Paired',guide = guide_legend(reverse = TRUE))+
# scale_fill_brewer(type = 'qual', palette = 'Paired',guide = guide_legend(reverse = TRUE))+
scale_color_viridis(discrete = TRUE, guide = guide_legend(reverse = TRUE)) +
scale_fill_viridis(discrete = TRUE, guide = guide_legend(reverse = TRUE)) +
scale_y_discrete(limits = rev) +
coord_fixed(ratio = 0.05, xlim = c(-0.23, 0.15)) +
# scale_shape(guide = guide_legend(reverse = TRUE)) +
# scale_x_continuous(limits = c(-2, 2), expand = c(0, 0), breaks = seq(-2,2, by = 1), labels = seq(-2,2, by = 1)) +
labs(y = NULL, x = "Google Mobility\n[Standardized effect sizes on\nflight initiation distances]", tag = "c)") + # title = "b) Effect of ") +
# ylim(c(0,100))+
# coord_flip()+
theme_bw() +
theme(
legend.position = "none",
plot.subtitle = element_text(size = 7),
plot.title = element_text(size = 7),
plot.tag = element_text(size = 7),
legend.title = element_text(size = 7),
legend.text = element_text(size = 6),
## legend.spacing.y = unit(0.1, 'cm'),
legend.key.height = unit(0.5, "line"),
# plot.margin = margin(b = 0.5, l = 0.5, t = 0.5, r =0.5, unit = "pt"),
panel.grid = element_blank(),
panel.border = element_blank(),
panel.background = element_blank(),
axis.line = element_line(colour = ax_lines, size = 0.25),
axis.line.y = element_blank(),
axis.ticks.y = element_blank(),
axis.ticks.x = element_line(colour = ax_lines, size = 0.25),
axis.ticks.length = unit(1, "pt"),
axis.text.x = element_text(colour = "black", size = 6),
axis.text.y = element_text(colour = "black", size = 7),
axis.title = element_text(size = 7)
)
#gg0
# ggsave(here::here('Outputs/Figure_Sz.png'),g0, width = 30, height =5, units = 'cm')
# combine
grid.draw(rbind(ggplotGrob(gs2_), ggplotGrob(g0_), ggplotGrob(gg0_)))
if(save_plot==TRUE){
ggsave(here::here("Outputs/Fig_S2_rev_v6.png"), rbind(ggplotGrob(gs2_), ggplotGrob(g0_), ggplotGrob(gg0_)), width = 30, height = 15, units = "cm")
}Figure S1 | Comparing estimates from alternative models. Changes in avian tolerance towards humans in response to (a) Period (before vs during the COVID-19 shutdowns) (b) stringency of governmental measures and (c) Google Mobility. The dots with horizontal lines represent the estimated standardised effect size and their 95% confidence intervals based on the joint posterior distribution of 5,000 simulated values generated by the sim function from the arm R-package (Gelman et al. 2016) from the output of the mixed models (for details see Table S2). The name of each effect size highlights the corresponding model in Table S2a for (a), Table S2b for (b) and Table S2c for (c), the random structure of the specific model, if applicable, the condition used to reduce the dataset, and sample size. Depicted are effect sizes based on full (01) and reduced datasets with ≥5 (02) or ≥10 observations per species and period (03). For plots of model assumptions see https://doi.org/10.17605/OSF.IO/WUZH7 (Bulla et al. 2022). Note that effect sizes are small and estimates centre around zero.
Table S2a | Alternative models on escape distance given Period
m1a_ = m_out(name = "Table S2a - 1a", dep = "Escape distance", model = m1a, nsim = 5000)
m1b_ = m_out(name = "Table S2a - 1b", dep = "Escape distance", model = m1b, nsim = 5000)
m1c_ = m_out(name = "Table S2a - 1c", dep = "Escape distance", model = m1c, nsim = 5000)
m1d_ = m_out(name = "Table S2a - 1d", dep = "Escape distance", model = m1d, nsim = 5000)
m1e_ = m_out(name = "Table S2a - 1e", dep = "Escape distance", model = m1e, nsim = 5000)
m2a_ = m_out(name = "Table S2a - 2a", dep = "Escape distance", model = m2a, nsim = 5000)
m2b_ = m_out(name = "Table S2a - 2b", dep = "Escape distance", model = m2b, nsim = 5000)
m2c_ = m_out(name = "Table S2a - 2c", dep = "Escape distance", model = m2c, nsim = 5000)
m3a_ = m_out(name = "Table S2a - 3a", dep = "Escape distance", model = m3a, nsim = 5000)
m3b_ = m_out(name = "Table S2a - 3b", dep = "Escape distance", model = m3b, nsim = 5000)
m3c_ = m_out(name = "Table S2a - 3c", dep = "Escape distance", model = m3c, nsim = 5000)
out1 = rbind(m1a_, m1b_, m1c_, m1d_, m1e_, m2a_, m2b_, m2c_, m3a_, m3b_, m3c_, fill = TRUE)
out1[is.na(out1)] = ""
out1[, effect := gsub("scale\\(Covid\\)", "Period", effect)]
out1[, effect := gsub("scale\\(Year\\)", "year", effect)]
out1[, effect := gsub("scale\\(log\\(SD\\)\\)", "starting distance (ln)", effect)]
out1[, effect := gsub("scale\\(Temp\\)", "temperaturre", effect)]
out1[, effect := gsub("scale\\(log\\(FlockSize\\)\\)", "flock size (ln)", effect)]
out1[, effect := gsub("scale\\(log\\(BodyMass\\)\\)", "body mass (ln)", effect)]
out1[, effect := gsub("scale\\(sin\\(rad\\)\\)", "time (sine of radians)", effect)]
out1[, effect := gsub("scale\\(cos\\(rad\\)\\)", "time (cosine of radians)", effect)]
out1[, effect := gsub("Species", "species", effect)]
out1[, effect := gsub("Country", "country", effect)]
out1[, effect := gsub("Year", "year", effect)]
out1[, effect := gsub("sp_day_year", "species within day & year", effect)]
out1[, effect := gsub("IDLocality", "site", effect)]
out1[, effect := gsub("sp_loc", "species within site", effect)]
out1[type == "random" & grepl("Period", effect, fixed = TRUE), effect := paste("Period (slope) |", gsub(" Period", "", effect))]
out1$R2_mar=out1$R2_con=NULL
fwrite(file = here::here("Outputs/Table_S2a.csv"), out1)
out1$response = out1$error_structure = NULL
out1[model != "", model := paste0('0',substring(model, 12))]
setnames(out1, old = c("estimate_r", "lwr_r", "upr_r"), new = c("estimate", "lower", "upper"))
out1 %>%
kbl() %>%
kable_paper("hover", full_width = F)| model | N | type | effect | estimate | lower | upper |
|---|---|---|---|---|---|---|
| 0 1a | 6369 | fixed | (Intercept) | 0.082 | -0.129 | 0.286 |
| fixed | starting distance (ln) | 0.477 | 0.453 | 0.501 | ||
| fixed | flock size (ln) | -0.008 | -0.027 | 0.011 | ||
| fixed | body mass (ln) | 0.048 | -0.02 | 0.111 | ||
| fixed | time (sine of radians) | 0.022 | -0.006 | 0.05 | ||
| fixed | time (cosine of radians) | 0.028 | 0.006 | 0.051 | ||
| fixed | temperaturre | 0 | -0.029 | 0.028 | ||
| fixed | Period | -0.033 | -0.198 | 0.132 | ||
| random | species within day & year (Intercept) | 8% | 6% | 9% | ||
| random | species within site (Intercept) | 5% | 4% | 5% | ||
| random | site (Intercept) | 1% | 0% | 5% | ||
| random | Period (slope) | site | 1% | 0% | 5% | ||
| random | species (Intercept) | 9% | 8% | 8% | ||
| random | genus (Intercept) | 6% | 6% | 6% | ||
| random | weekday (Intercept) | 0% | 0% | 0% | ||
| random | country (Intercept) | 1% | -3% | 5% | ||
| random | Period (slope) | country | 1% | -3% | 5% | ||
| random | year (Intercept) | 3% | 1% | 4% | ||
| random | Residual | 65% | 51% | 76% | ||
| 0 1b | 6369 | fixed | (Intercept) | 0.068 | -0.155 | 0.294 |
| fixed | starting distance (ln) | 0.478 | 0.455 | 0.502 | ||
| fixed | flock size (ln) | -0.007 | -0.026 | 0.012 | ||
| fixed | body mass (ln) | 0.048 | -0.017 | 0.113 | ||
| fixed | time (sine of radians) | 0.021 | -0.006 | 0.049 | ||
| fixed | time (cosine of radians) | 0.027 | 0.004 | 0.05 | ||
| fixed | temperaturre | -0.002 | -0.03 | 0.028 | ||
| fixed | Period | -0.048 | -0.204 | 0.113 | ||
| random | species within day & year (Intercept) | 12% | 7% | 14% | ||
| random | species within site (Intercept) | 0% | -1% | 0% | ||
| random | Period (slope) | species within site | 0% | -1% | 0% | ||
| random | site (Intercept) | 2% | 0% | 5% | ||
| random | Period (slope) | site | 2% | 0% | 5% | ||
| random | species (Intercept) | 0% | -1% | 8% | ||
| random | Period (slope) | species | 0% | -1% | 8% | ||
| random | genus (Intercept) | 0% | 0% | 5% | ||
| random | Period (slope) | genus | 0% | 0% | 5% | ||
| random | weekday (Intercept) | 0% | 0% | 0% | ||
| random | country (Intercept) | 1% | -4% | 6% | ||
| random | Period (slope) | country | 1% | -4% | 6% | ||
| random | year (Intercept) | 3% | 2% | 4% | ||
| random | Residual | 79% | 42% | 95% | ||
| 0 1c | 6369 | fixed | (Intercept) | 0.015 | -0.148 | 0.181 |
| fixed | starting distance (ln) | 0.494 | 0.469 | 0.517 | ||
| fixed | flock size (ln) | -0.018 | -0.037 | 0 | ||
| fixed | body mass (ln) | 0.024 | -0.007 | 0.054 | ||
| fixed | time (sine of radians) | 0.015 | -0.014 | 0.043 | ||
| fixed | time (cosine of radians) | 0.027 | 0.004 | 0.051 | ||
| fixed | temperaturre | -0.004 | -0.034 | 0.026 | ||
| fixed | Period | -0.041 | -0.212 | 0.136 | ||
| random | species within day & year (Intercept) | 12% | 7% | 13% | ||
| random | species within site (Intercept) | 1% | 0% | 12% | ||
| random | Period (slope) | species within site | 1% | 0% | 12% | ||
| random | site (Intercept) | 1% | 1% | 6% | ||
| random | Period (slope) | site | 1% | 1% | 6% | ||
| random | Period (slope) | species | 1% | 1% | 1% | ||
| random | Period (slope) | genus | 0% | 0% | 0% | ||
| random | weekday (Intercept) | 0% | 0% | 0% | ||
| random | country (Intercept) | 1% | -2% | 4% | ||
| random | Period (slope) | country | 1% | -2% | 4% | ||
| random | year (Intercept) | 3% | 2% | 4% | ||
| random | Residual | 77% | 46% | 87% | ||
| 0 1d | 6369 | fixed | (Intercept) | 0.084 | -0.128 | 0.292 |
| fixed | starting distance (ln) | 0.476 | 0.452 | 0.5 | ||
| fixed | flock size (ln) | -0.008 | -0.027 | 0.011 | ||
| fixed | body mass (ln) | 0.048 | -0.017 | 0.116 | ||
| fixed | time (sine of radians) | 0.022 | -0.006 | 0.05 | ||
| fixed | time (cosine of radians) | 0.028 | 0.004 | 0.051 | ||
| fixed | temperaturre | 0 | -0.029 | 0.029 | ||
| fixed | Period | -0.035 | -0.2 | 0.132 | ||
| random | species within day & year (Intercept) | 8% | 6% | 9% | ||
| random | species within site (Intercept) | 0% | 0% | 4% | ||
| random | Period (slope) | species within site | 0% | 0% | 4% | ||
| random | site (Intercept) | 1% | 0% | 5% | ||
| random | Period (slope) | site | 1% | 0% | 5% | ||
| random | species (Intercept) | 10% | 9% | 10% | ||
| random | genus (Intercept) | 7% | 6% | 6% | ||
| random | weekday (Intercept) | 0% | 0% | 0% | ||
| random | country (Intercept) | 1% | -3% | 5% | ||
| random | Period (slope) | country | 1% | -3% | 5% | ||
| random | year (Intercept) | 3% | 2% | 4% | ||
| random | Residual | 66% | 48% | 78% | ||
| 0 1e | 6369 | fixed | (Intercept) | 0.075 | -0.125 | 0.282 |
| fixed | starting distance (ln) | 0.478 | 0.454 | 0.5 | ||
| fixed | flock size (ln) | -0.008 | -0.027 | 0.011 | ||
| fixed | body mass (ln) | 0.048 | -0.019 | 0.116 | ||
| fixed | time (sine of radians) | 0.022 | -0.007 | 0.049 | ||
| fixed | time (cosine of radians) | 0.024 | 0.001 | 0.047 | ||
| fixed | temperaturre | -0.006 | -0.034 | 0.023 | ||
| fixed | Period | -0.033 | -0.203 | 0.135 | ||
| random | species within day & year (Intercept) | 7% | 6% | 8% | ||
| random | species within site (Intercept) | 5% | 4% | 5% | ||
| random | site (Intercept) | 6% | 6% | 6% | ||
| random | species (Intercept) | 8% | 8% | 9% | ||
| random | genus (Intercept) | 6% | 5% | 6% | ||
| random | weekday (Intercept) | 0% | 0% | 1% | ||
| random | country (Intercept) | 2% | -3% | 5% | ||
| random | Period (slope) | country | 2% | -3% | 5% | ||
| random | year (Intercept) | 2% | 1% | 4% | ||
| random | Residual | 62% | 52% | 71% | ||
| 0 2a | 5261 | fixed | (Intercept) | 0.159 | -0.183 | 0.512 |
| fixed | starting distance (ln) | 0.462 | 0.437 | 0.489 | ||
| fixed | flock size (ln) | -0.01 | -0.031 | 0.011 | ||
| fixed | body mass (ln) | -0.014 | -0.124 | 0.096 | ||
| fixed | time (sine of radians) | 0.028 | -0.001 | 0.057 | ||
| fixed | time (cosine of radians) | 0.036 | 0.012 | 0.06 | ||
| fixed | temperaturre | -0.006 | -0.038 | 0.027 | ||
| fixed | Period | -0.036 | -0.211 | 0.135 | ||
| random | species within day & year (Intercept) | 7% | 4% | 9% | ||
| random | species within site (Intercept) | 3% | 2% | 4% | ||
| random | site (Intercept) | 1% | -1% | 2% | ||
| random | Period (slope) | site | 1% | -1% | 2% | ||
| random | species (Intercept) | 10% | 8% | 9% | ||
| random | genus (Intercept) | 5% | 4% | 4% | ||
| random | weekday (Intercept) | 0% | 0% | 0% | ||
| random | country (Intercept) | 2% | -6% | 20% | ||
| random | Period (slope) | country | 2% | -6% | 20% | ||
| random | year (Intercept) | 3% | 2% | 3% | ||
| random | Residual | 64% | 34% | 84% | ||
| 0 2b | 5261 | fixed | (Intercept) | 0.16 | -0.192 | 0.493 |
| fixed | starting distance (ln) | 0.462 | 0.436 | 0.489 | ||
| fixed | flock size (ln) | -0.01 | -0.031 | 0.01 | ||
| fixed | body mass (ln) | -0.011 | -0.12 | 0.1 | ||
| fixed | time (sine of radians) | 0.027 | -0.002 | 0.056 | ||
| fixed | time (cosine of radians) | 0.036 | 0.011 | 0.061 | ||
| fixed | temperaturre | -0.005 | -0.037 | 0.026 | ||
| fixed | Period | -0.035 | -0.211 | 0.137 | ||
| random | species within day & year (Intercept) | 9% | 4% | 11% | ||
| random | species within site (Intercept) | 0% | 0% | 2% | ||
| random | Period (slope) | species within site | 0% | 0% | 2% | ||
| random | site (Intercept) | 1% | -1% | 2% | ||
| random | Period (slope) | site | 1% | -1% | 2% | ||
| random | species (Intercept) | 0% | -1% | 5% | ||
| random | Period (slope) | species | 0% | -1% | 5% | ||
| random | genus (Intercept) | 0% | 0% | 6% | ||
| random | Period (slope) | genus | 0% | 0% | 6% | ||
| random | weekday (Intercept) | 0% | 0% | 0% | ||
| random | country (Intercept) | 3% | -6% | 17% | ||
| random | Period (slope) | country | 3% | -6% | 17% | ||
| random | year (Intercept) | 4% | 2% | 3% | ||
| random | Residual | 78% | 30% | 100% | ||
| 0 2c | 5261 | fixed | (Intercept) | 0.152 | -0.183 | 0.503 |
| fixed | starting distance (ln) | 0.464 | 0.438 | 0.491 | ||
| fixed | flock size (ln) | -0.01 | -0.03 | 0.01 | ||
| fixed | body mass (ln) | -0.008 | -0.121 | 0.099 | ||
| fixed | time (sine of radians) | 0.026 | -0.002 | 0.056 | ||
| fixed | time (cosine of radians) | 0.032 | 0.008 | 0.057 | ||
| fixed | temperaturre | -0.012 | -0.043 | 0.02 | ||
| fixed | Period | -0.032 | -0.217 | 0.145 | ||
| random | species within day & year (Intercept) | 8% | 5% | 9% | ||
| random | species within site (Intercept) | 3% | 2% | 4% | ||
| random | site (Intercept) | 4% | 3% | 4% | ||
| random | species (Intercept) | 8% | 7% | 7% | ||
| random | genus (Intercept) | 6% | 4% | 6% | ||
| random | weekday (Intercept) | 0% | 0% | 0% | ||
| random | country (Intercept) | 2% | -6% | 17% | ||
| random | Period (slope) | country | 2% | -6% | 17% | ||
| random | year (Intercept) | 3% | 2% | 4% | ||
| random | Residual | 63% | 39% | 81% | ||
| 0 3a | 5107 | fixed | (Intercept) | 0.133 | -0.216 | 0.476 |
| fixed | starting distance (ln) | 0.462 | 0.435 | 0.489 | ||
| fixed | flock size (ln) | -0.01 | -0.03 | 0.011 | ||
| fixed | body mass (ln) | 0.006 | -0.104 | 0.113 | ||
| fixed | time (sine of radians) | 0.027 | -0.003 | 0.057 | ||
| fixed | time (cosine of radians) | 0.034 | 0.009 | 0.06 | ||
| fixed | temperaturre | -0.007 | -0.04 | 0.026 | ||
| fixed | Period | -0.028 | -0.196 | 0.144 | ||
| random | species within day & year (Intercept) | 7% | 4% | 9% | ||
| random | species within site (Intercept) | 3% | 2% | 4% | ||
| random | site (Intercept) | 1% | -1% | 2% | ||
| random | Period (slope) | site | 1% | -1% | 2% | ||
| random | species (Intercept) | 21% | 16% | 18% | ||
| random | genus (Intercept) | 0% | 0% | 0% | ||
| random | weekday (Intercept) | 0% | 0% | 0% | ||
| random | country (Intercept) | 2% | -5% | 19% | ||
| random | Period (slope) | country | 2% | -5% | 19% | ||
| random | year (Intercept) | 3% | 2% | 3% | ||
| random | Residual | 60% | 32% | 79% | ||
| 0 3b | 5107 | fixed | (Intercept) | 0.141 | -0.197 | 0.491 |
| fixed | starting distance (ln) | 0.461 | 0.433 | 0.489 | ||
| fixed | flock size (ln) | -0.01 | -0.03 | 0.011 | ||
| fixed | body mass (ln) | 0.013 | -0.104 | 0.13 | ||
| fixed | time (sine of radians) | 0.027 | -0.003 | 0.058 | ||
| fixed | time (cosine of radians) | 0.034 | 0.009 | 0.059 | ||
| fixed | temperaturre | -0.006 | -0.04 | 0.026 | ||
| fixed | Period | -0.03 | -0.203 | 0.141 | ||
| random | species within day & year (Intercept) | 9% | 3% | 11% | ||
| random | species within site (Intercept) | 0% | 0% | 2% | ||
| random | Period (slope) | species within site | 0% | 0% | 2% | ||
| random | site (Intercept) | 1% | -1% | 2% | ||
| random | Period (slope) | site | 1% | -1% | 2% | ||
| random | species (Intercept) | 0% | -1% | 9% | ||
| random | Period (slope) | species | 0% | -1% | 9% | ||
| random | genus (Intercept) | 0% | 0% | 3% | ||
| random | Period (slope) | genus | 0% | 0% | 3% | ||
| random | weekday (Intercept) | 0% | 0% | 0% | ||
| random | country (Intercept) | 3% | -5% | 16% | ||
| random | Period (slope) | country | 3% | -5% | 16% | ||
| random | year (Intercept) | 4% | 2% | 3% | ||
| random | Residual | 79% | 28% | 100% | ||
| 0 3c | 5107 | fixed | (Intercept) | 0.127 | -0.221 | 0.478 |
| fixed | starting distance (ln) | 0.464 | 0.437 | 0.49 | ||
| fixed | flock size (ln) | -0.01 | -0.031 | 0.01 | ||
| fixed | body mass (ln) | 0.005 | -0.105 | 0.12 | ||
| fixed | time (sine of radians) | 0.027 | -0.003 | 0.056 | ||
| fixed | time (cosine of radians) | 0.028 | 0.004 | 0.053 | ||
| fixed | temperaturre | -0.014 | -0.046 | 0.018 | ||
| fixed | Period | -0.031 | -0.208 | 0.145 | ||
| random | species within day & year (Intercept) | 7% | 5% | 9% | ||
| random | species within site (Intercept) | 3% | 2% | 4% | ||
| random | site (Intercept) | 5% | 3% | 5% | ||
| random | species (Intercept) | 15% | 11% | 15% | ||
| random | genus (Intercept) | 0% | 0% | 0% | ||
| random | weekday (Intercept) | 0% | 0% | 0% | ||
| random | country (Intercept) | 2% | -6% | 18% | ||
| random | Period (slope) | country | 2% | -6% | 18% | ||
| random | year (Intercept) | 3% | 2% | 4% | ||
| random | Residual | 62% | 39% | 78% |
Note that (1a) model is the one reported in the main text.
Table S2b | Alternative models on escape distance given Stringency
m01a_ = m_out(name = "Table S2b - 1a", dep = "Escape distance", model = m01a, nsim = 5000)
m01b_ = m_out(name = "Table S2b - 1b", dep = "Escape distance", model = m01b, nsim = 5000)
m01c_ = m_out(name = "Table S2b - 1c", dep = "Escape distance", model = m01c, nsim = 5000)
m02a_ = m_out(name = "Table S2b - 2a", dep = "Escape distance", model = m02a, nsim = 5000)
m02b_ = m_out(name = "Table S2b - 2b", dep = "Escape distance", model = m02b, nsim = 5000)
m02c_ = m_out(name = "Table S2b - 2c", dep = "Escape distance", model = m02c, nsim = 5000)
m03a_ = m_out(name = "Table S2b - 3a", dep = "Escape distance", model = m03a, nsim = 5000)
m03b_ = m_out(name = "Table S2b - 3b", dep = "Escape distance", model = m03b, nsim = 5000)
m03c_ = m_out(name = "Table S2b - 3c", dep = "Escape distance", model = m03c, nsim = 5000)
out2 = rbind(m01a_, m01b_, m01c_, m02a_, m02b_, m02c_, m03a_, m03b_, m03c_, fill = TRUE)
out2[is.na(out2)] = ""
out2[, effect := gsub("scale\\(StringencyIndex\\)", "stringency index", effect)]
out2[, effect := gsub("scale\\(Year\\)", "year", effect)]
out2[, effect := gsub("scale\\(log\\(SD\\)\\)", "starting distance (ln)", effect)]
out2[, effect := gsub("scale\\(Temp\\)", "temperaturre", effect)]
out2[, effect := gsub("scale\\(log\\(FlockSize\\)\\)", "flock size (ln)", effect)]
out2[, effect := gsub("scale\\(log\\(BodyMass\\)\\)", "body mass (ln)", effect)]
out2[, effect := gsub("scale\\(sin\\(rad\\)\\)", "time (sine of radians)", effect)]
out2[, effect := gsub("scale\\(cos\\(rad\\)\\)", "time (cosine of radians)", effect)]
out2[, effect := gsub("Species", "species", effect)]
out2[, effect := gsub("Country", "country", effect)]
out2[, effect := gsub("Year", "year", effect)]
out2[, effect := gsub("sp_day_year", "species within day & year", effect)]
out2[, effect := gsub("IDLocality", "site", effect)]
out2[, effect := gsub("sp_loc", "species within site", effect)]
out2[type == "random" & grepl("stringency index", effect, fixed = TRUE), effect := paste("stringency index (slope) |", gsub(" stringency index", "", effect))]
out2$R2_mar = out2$R2_con = NULL
fwrite(file = here::here("Outputs/Table_S2b.csv"), out2)
out2$response = out2$error_structure = NULL
out2[model != "", model := paste0("0", substring(model, 12))]
setnames(out2, old = c("estimate_r", "lwr_r", "upr_r"), new = c("estimate", "lower", "upper"))
out2 %>%
kbl() %>%
kable_paper("hover", full_width = F)| model | N | type | effect | estimate | lower | upper |
|---|---|---|---|---|---|---|
| 0 1a | 3676 | fixed | (Intercept) | -0.074 | -0.345 | 0.173 |
| fixed | year | 0.02 | -0.028 | 0.068 | ||
| fixed | starting distance (ln) | 0.492 | 0.46 | 0.524 | ||
| fixed | flock size (ln) | -0.003 | -0.03 | 0.023 | ||
| fixed | body mass (ln) | 0.021 | -0.046 | 0.091 | ||
| fixed | time (sine of radians) | -0.02 | -0.061 | 0.022 | ||
| fixed | time (cosine of radians) | -0.005 | -0.04 | 0.032 | ||
| fixed | temperaturre | -0.007 | -0.049 | 0.036 | ||
| fixed | stringency index | 0.018 | -0.157 | 0.182 | ||
| random | species within day & year (Intercept) | 5% | 3% | 7% | ||
| random | species within site (Intercept) | 5% | 3% | 6% | ||
| random | species (Intercept) | 7% | 6% | 9% | ||
| random | site (Intercept) | 8% | 7% | 9% | ||
| random | genus (Intercept) | 4% | 4% | 4% | ||
| random | weekday (Intercept) | 0% | 0% | 0% | ||
| random | country (Intercept) | 2% | -16% | 16% | ||
| random | stringency index (slope) | country | 2% | -16% | 16% | ||
| random | Residual | 67% | 45% | 96% | ||
| 0 1b | 3676 | fixed | (Intercept) | -0.073 | -0.319 | 0.173 |
| fixed | year | 0.019 | -0.028 | 0.065 | ||
| fixed | starting distance (ln) | 0.493 | 0.462 | 0.525 | ||
| fixed | flock size (ln) | -0.003 | -0.029 | 0.023 | ||
| fixed | body mass (ln) | 0.017 | -0.052 | 0.085 | ||
| fixed | time (sine of radians) | -0.02 | -0.061 | 0.021 | ||
| fixed | time (cosine of radians) | -0.004 | -0.039 | 0.03 | ||
| fixed | temperaturre | -0.008 | -0.049 | 0.033 | ||
| fixed | stringency index | 0.006 | -0.16 | 0.174 | ||
| random | species within day & year (Intercept) | 5% | 3% | 7% | ||
| random | species within site (Intercept) | 5% | 3% | 7% | ||
| random | species (Intercept) | 7% | 5% | 8% | ||
| random | site (Intercept) | 8% | 6% | 10% | ||
| random | genus (Intercept) | 0% | -1% | 4% | ||
| random | stringency index (slope) | genus | 0% | -1% | 4% | ||
| random | weekday (Intercept) | 0% | 0% | 0% | ||
| random | country (Intercept) | 2% | -16% | 15% | ||
| random | stringency index (slope) | country | 2% | -16% | 15% | ||
| random | Residual | 70% | 44% | 100% | ||
| 0 1c | 3676 | fixed | (Intercept) | -0.103 | -0.278 | 0.066 |
| fixed | year | -0.003 | -0.05 | 0.044 | ||
| fixed | starting distance (ln) | 0.549 | 0.518 | 0.581 | ||
| fixed | flock size (ln) | -0.023 | -0.05 | 0.003 | ||
| fixed | body mass (ln) | -0.024 | -0.053 | 0.005 | ||
| fixed | time (sine of radians) | -0.039 | -0.081 | 0.002 | ||
| fixed | time (cosine of radians) | 0.005 | -0.031 | 0.043 | ||
| fixed | temperaturre | 0.003 | -0.039 | 0.045 | ||
| fixed | stringency index | 0.021 | -0.14 | 0.181 | ||
| random | site (Intercept) | 10% | 10% | 10% | ||
| random | country (Intercept) | 1% | -7% | 8% | ||
| random | stringency index (slope) | country | 1% | -7% | 8% | ||
| random | Residual | 88% | 74% | 103% | ||
| 0 2a | 3573 | fixed | (Intercept) | -0.154 | -0.455 | 0.133 |
| fixed | year | 0.022 | -0.025 | 0.071 | ||
| fixed | starting distance (ln) | 0.483 | 0.451 | 0.516 | ||
| fixed | flock size (ln) | 0.003 | -0.023 | 0.029 | ||
| fixed | body mass (ln) | -0.024 | -0.1 | 0.052 | ||
| fixed | time (sine of radians) | -0.02 | -0.061 | 0.022 | ||
| fixed | time (cosine of radians) | -0.008 | -0.043 | 0.029 | ||
| fixed | temperaturre | -0.008 | -0.05 | 0.034 | ||
| fixed | stringency index | 0.024 | -0.153 | 0.191 | ||
| random | species within day & year (Intercept) | 5% | 3% | 7% | ||
| random | species within site (Intercept) | 4% | 3% | 7% | ||
| random | site (Intercept) | 7% | 5% | 9% | ||
| random | species (Intercept) | 11% | 8% | 13% | ||
| random | genus (Intercept) | 2% | 2% | 2% | ||
| random | weekday (Intercept) | 0% | 0% | 0% | ||
| random | country (Intercept) | 3% | -20% | 20% | ||
| random | stringency index (slope) | country | 3% | -20% | 20% | ||
| random | Residual | 65% | 39% | 103% | ||
| 0 2b | 3573 | fixed | (Intercept) | -0.151 | -0.441 | 0.146 |
| fixed | year | 0.021 | -0.029 | 0.07 | ||
| fixed | starting distance (ln) | 0.484 | 0.452 | 0.517 | ||
| fixed | flock size (ln) | 0.003 | -0.023 | 0.03 | ||
| fixed | body mass (ln) | -0.03 | -0.107 | 0.044 | ||
| fixed | time (sine of radians) | -0.02 | -0.06 | 0.021 | ||
| fixed | time (cosine of radians) | -0.007 | -0.043 | 0.03 | ||
| fixed | temperaturre | -0.008 | -0.051 | 0.034 | ||
| fixed | stringency index | 0.015 | -0.159 | 0.189 | ||
| random | species within day & year (Intercept) | 5% | 3% | 8% | ||
| random | species within site (Intercept) | 4% | 3% | 7% | ||
| random | site (Intercept) | 7% | 5% | 9% | ||
| random | species (Intercept) | 11% | 8% | 13% | ||
| random | genus (Intercept) | 0% | -1% | 2% | ||
| random | stringency index (slope) | genus | 0% | -1% | 2% | ||
| random | weekday (Intercept) | 0% | 0% | 0% | ||
| random | country (Intercept) | 3% | -21% | 19% | ||
| random | stringency index (slope) | country | 3% | -21% | 19% | ||
| random | Residual | 67% | 39% | 106% | ||
| 0 2c | 3573 | fixed | (Intercept) | -0.109 | -0.286 | 0.064 |
| fixed | year | -0.001 | -0.049 | 0.048 | ||
| fixed | starting distance (ln) | 0.543 | 0.511 | 0.576 | ||
| fixed | flock size (ln) | -0.018 | -0.045 | 0.009 | ||
| fixed | body mass (ln) | -0.027 | -0.055 | 0.003 | ||
| fixed | time (sine of radians) | -0.038 | -0.08 | 0.003 | ||
| fixed | time (cosine of radians) | 0.002 | -0.036 | 0.039 | ||
| fixed | temperaturre | 0.002 | -0.04 | 0.047 | ||
| fixed | stringency index | 0.023 | -0.136 | 0.19 | ||
| random | site (Intercept) | 10% | 10% | 10% | ||
| random | country (Intercept) | 1% | -8% | 8% | ||
| random | stringency index (slope) | country | 1% | -8% | 8% | ||
| random | Residual | 88% | 74% | 105% | ||
| 0 3a | 3425 | fixed | (Intercept) | -0.137 | -0.459 | 0.173 |
| fixed | year | 0.026 | -0.023 | 0.077 | ||
| fixed | starting distance (ln) | 0.48 | 0.448 | 0.512 | ||
| fixed | flock size (ln) | 0.001 | -0.027 | 0.028 | ||
| fixed | body mass (ln) | 0.023 | -0.061 | 0.11 | ||
| fixed | time (sine of radians) | -0.012 | -0.054 | 0.03 | ||
| fixed | time (cosine of radians) | -0.009 | -0.046 | 0.027 | ||
| fixed | temperaturre | -0.011 | -0.054 | 0.032 | ||
| fixed | stringency index | 0.02 | -0.153 | 0.187 | ||
| random | species within day & year (Intercept) | 5% | 3% | 7% | ||
| random | species within site (Intercept) | 5% | 3% | 7% | ||
| random | site (Intercept) | 8% | 5% | 10% | ||
| random | species (Intercept) | 9% | 7% | 10% | ||
| random | genus (Intercept) | 3% | 3% | 3% | ||
| random | weekday (Intercept) | 0% | 0% | 0% | ||
| random | country (Intercept) | 2% | -21% | 20% | ||
| random | stringency index (slope) | country | 2% | -21% | 20% | ||
| random | Residual | 67% | 39% | 105% | ||
| 0 3b | 3425 | fixed | (Intercept) | -0.14 | -0.456 | 0.181 |
| fixed | year | 0.027 | -0.022 | 0.075 | ||
| fixed | starting distance (ln) | 0.482 | 0.449 | 0.514 | ||
| fixed | flock size (ln) | 0.001 | -0.026 | 0.029 | ||
| fixed | body mass (ln) | 0.019 | -0.067 | 0.106 | ||
| fixed | time (sine of radians) | -0.011 | -0.053 | 0.032 | ||
| fixed | time (cosine of radians) | -0.009 | -0.045 | 0.028 | ||
| fixed | temperaturre | -0.012 | -0.055 | 0.033 | ||
| fixed | stringency index | 0.013 | -0.155 | 0.189 | ||
| random | species within day & year (Intercept) | 5% | 3% | 7% | ||
| random | species within site (Intercept) | 4% | 3% | 7% | ||
| random | site (Intercept) | 8% | 5% | 11% | ||
| random | species (Intercept) | 10% | 7% | 11% | ||
| random | genus (Intercept) | 0% | 0% | 2% | ||
| random | stringency index (slope) | genus | 0% | 0% | 2% | ||
| random | weekday (Intercept) | 0% | 0% | 0% | ||
| random | country (Intercept) | 3% | -23% | 21% | ||
| random | stringency index (slope) | country | 3% | -23% | 21% | ||
| random | Residual | 67% | 36% | 110% | ||
| 0 3c | 3425 | fixed | (Intercept) | -0.107 | -0.288 | 0.076 |
| fixed | year | 0.005 | -0.045 | 0.054 | ||
| fixed | starting distance (ln) | 0.541 | 0.508 | 0.574 | ||
| fixed | flock size (ln) | -0.02 | -0.047 | 0.007 | ||
| fixed | body mass (ln) | -0.02 | -0.048 | 0.01 | ||
| fixed | time (sine of radians) | -0.03 | -0.073 | 0.012 | ||
| fixed | time (cosine of radians) | -0.001 | -0.038 | 0.037 | ||
| fixed | temperaturre | -0.002 | -0.045 | 0.042 | ||
| fixed | stringency index | 0.023 | -0.137 | 0.183 | ||
| random | site (Intercept) | 11% | 11% | 11% | ||
| random | country (Intercept) | 1% | -8% | 8% | ||
| random | stringency index (slope) | country | 1% | -8% | 8% | ||
| random | Residual | 87% | 73% | 105% |
Note that (1a) model is the one reported in the main text.
Table S2c | Alternative models on escape distance given Google Mobility
# google
mg01a_ = m_out(name = "Table S2c - 1a", dep = "Escape distance", model = mg01a, nsim = 5000)
mg01b_ = m_out(name = "Table S2c - 1b", dep = "Escape distance", model = mg01b, nsim = 5000)
mg01c_ = m_out(name = "Table S2c - 1c", dep = "Escape distance", model = mg01c, nsim = 5000)
mg02a_ = m_out(name = "Table S2c - 2a", dep = "Escape distance", model = mg02a, nsim = 5000)
mg02b_ = m_out(name = "Table S2c - 2b", dep = "Escape distance", model = mg02b, nsim = 5000)
mg02c_ = m_out(name = "Table S2c - 2c", dep = "Escape distance", model = mg02c, nsim = 5000)
mg03a_ = m_out(name = "Table S2c - 3a", dep = "Escape distance", model = mg03a, nsim = 5000)
mg03b_ = m_out(name = "Table S2c - 3b", dep = "Escape distance", model = mg03b, nsim = 5000)
mg03c_ = m_out(name = "Table S2c - 3c", dep = "Escape distancey", model = mg03c, nsim = 5000)
out3 = rbind(mg01a_, mg01b_, mg01c_, mg02a_, mg02b_, mg02c_, mg03a_, mg03b_, mg03c_, fill = TRUE)
out3[is.na(out3)] = ""
out3[, effect := gsub("scale\\(parks_percent_change_from_baseline\\)", "Google Mobility", effect)]
out3[, effect := gsub("scale\\(Year\\)", "year", effect)]
out3[, effect := gsub("scale\\(log\\(SD\\)\\)", "starting distance (ln)", effect)]
out3[, effect := gsub("scale\\(Temp\\)", "temperaturre", effect)]
out3[, effect := gsub("scale\\(log\\(FlockSize\\)\\)", "flock size (ln)", effect)]
out3[, effect := gsub("scale\\(log\\(BodyMass\\)\\)", "body mass (ln)", effect)]
out3[, effect := gsub("scale\\(sin\\(rad\\)\\)", "time (sine of radians)", effect)]
out3[, effect := gsub("scale\\(cos\\(rad\\)\\)", "time (cosine of radians)", effect)]
out3[, effect := gsub("Species", "species", effect)]
out3[, effect := gsub("Country", "country", effect)]
out3[, effect := gsub("Year", "year", effect)]
out3[, effect := gsub("sp_day_year", "species within day & year", effect)]
out3[, effect := gsub("IDLocality", "site", effect)]
out3[, effect := gsub("sp_loc", "species within site", effect)]
out3[type == "random" & grepl("Google Mobility", effect, fixed = TRUE), effect := paste("Google Mobility (slope) |", gsub(" Google Mobility", "", effect))]
out3$R2_mar = out3$R2_con = NULL
fwrite(file = here::here("Outputs/Table_S2c.csv"), out3)
out3$response = out3$error_structure = NULL
out3[model != "", model := paste0("0", substring(model, 12))]
setnames(out3, old = c("estimate_r", "lwr_r", "upr_r"), new = c("estimate", "lower", "upper"))
out3 %>%
kbl() %>%
kable_paper("hover", full_width = F)| model | N | type | effect | estimate | lower | upper |
|---|---|---|---|---|---|---|
| 0 1a | 3644 | fixed | (Intercept) | 0.025 | -0.097 | 0.146 |
| fixed | year | 0.019 | -0.023 | 0.06 | ||
| fixed | starting distance (ln) | 0.499 | 0.468 | 0.531 | ||
| fixed | flock size (ln) | -0.001 | -0.028 | 0.026 | ||
| fixed | body mass (ln) | 0.025 | -0.042 | 0.092 | ||
| fixed | time (sine of radians) | -0.009 | -0.05 | 0.033 | ||
| fixed | time (cosine of radians) | 0.01 | -0.023 | 0.043 | ||
| fixed | temperaturre | 0.021 | -0.017 | 0.059 | ||
| fixed | Google Mobility | -0.1 | -0.227 | 0.026 | ||
| random | species within day & year (Intercept) | 6% | 5% | 6% | ||
| random | species within site (Intercept) | 5% | 4% | 5% | ||
| random | species (Intercept) | 7% | 7% | 8% | ||
| random | site (Intercept) | 9% | 8% | 10% | ||
| random | genus (Intercept) | 3% | 3% | 4% | ||
| random | weekday (Intercept) | 0% | 0% | 0% | ||
| random | country (Intercept) | 1% | -3% | 4% | ||
| random | Google Mobility (slope) | country | 1% | -3% | 4% | ||
| random | Residual | 68% | 61% | 78% | ||
| 0 1b | 3644 | fixed | (Intercept) | 0.024 | -0.093 | 0.141 |
| fixed | year | 0.019 | -0.021 | 0.061 | ||
| fixed | starting distance (ln) | 0.498 | 0.466 | 0.53 | ||
| fixed | flock size (ln) | -0.001 | -0.027 | 0.025 | ||
| fixed | body mass (ln) | 0.029 | -0.038 | 0.095 | ||
| fixed | time (sine of radians) | -0.01 | -0.05 | 0.032 | ||
| fixed | time (cosine of radians) | 0.01 | -0.023 | 0.042 | ||
| fixed | temperaturre | 0.022 | -0.017 | 0.06 | ||
| fixed | Google Mobility | -0.101 | -0.232 | 0.026 | ||
| random | species within day & year (Intercept) | 6% | 5% | 6% | ||
| random | species within site (Intercept) | 5% | 4% | 5% | ||
| random | species (Intercept) | 7% | 7% | 7% | ||
| random | site (Intercept) | 9% | 9% | 9% | ||
| random | genus (Intercept) | 0% | -1% | 4% | ||
| random | Google Mobility (slope) | genus | 0% | -1% | 4% | ||
| random | weekday (Intercept) | 0% | 0% | 0% | ||
| random | country (Intercept) | 1% | -3% | 4% | ||
| random | Google Mobility (slope) | country | 1% | -3% | 4% | ||
| random | Residual | 71% | 59% | 81% | ||
| 0 1c | 3644 | fixed | (Intercept) | -0.057 | -0.143 | 0.026 |
| fixed | year | 0.013 | -0.027 | 0.053 | ||
| fixed | starting distance (ln) | 0.559 | 0.528 | 0.59 | ||
| fixed | flock size (ln) | -0.021 | -0.048 | 0.004 | ||
| fixed | body mass (ln) | -0.023 | -0.052 | 0.005 | ||
| fixed | time (sine of radians) | -0.033 | -0.075 | 0.009 | ||
| fixed | time (cosine of radians) | 0.02 | -0.014 | 0.056 | ||
| fixed | temperaturre | 0.031 | -0.007 | 0.069 | ||
| fixed | Google Mobility | -0.087 | -0.206 | 0.034 | ||
| random | site (Intercept) | 11% | 10% | 12% | ||
| random | country (Intercept) | 0% | -1% | 2% | ||
| random | Google Mobility (slope) | country | 0% | -1% | 2% | ||
| random | Residual | 88% | 83% | 92% | ||
| 0 2a | 3545 | fixed | (Intercept) | -0.027 | -0.171 | 0.11 |
| fixed | year | 0.019 | -0.025 | 0.061 | ||
| fixed | starting distance (ln) | 0.491 | 0.458 | 0.522 | ||
| fixed | flock size (ln) | 0.004 | -0.023 | 0.031 | ||
| fixed | body mass (ln) | -0.019 | -0.094 | 0.054 | ||
| fixed | time (sine of radians) | -0.01 | -0.051 | 0.032 | ||
| fixed | time (cosine of radians) | 0.009 | -0.025 | 0.043 | ||
| fixed | temperaturre | 0.019 | -0.018 | 0.059 | ||
| fixed | Google Mobility | -0.098 | -0.232 | 0.034 | ||
| random | species within day & year (Intercept) | 6% | 5% | 6% | ||
| random | species within site (Intercept) | 4% | 4% | 5% | ||
| random | site (Intercept) | 9% | 8% | 9% | ||
| random | species (Intercept) | 10% | 9% | 11% | ||
| random | genus (Intercept) | 2% | 1% | 2% | ||
| random | weekday (Intercept) | 0% | 0% | 0% | ||
| random | country (Intercept) | 0% | -5% | 5% | ||
| random | Google Mobility (slope) | country | 0% | -5% | 5% | ||
| random | Residual | 69% | 59% | 80% | ||
| 0 2b | 3545 | fixed | (Intercept) | -0.035 | -0.168 | 0.098 |
| fixed | year | 0.019 | -0.024 | 0.061 | ||
| fixed | starting distance (ln) | 0.49 | 0.458 | 0.522 | ||
| fixed | flock size (ln) | 0.004 | -0.022 | 0.03 | ||
| fixed | body mass (ln) | -0.03 | -0.101 | 0.043 | ||
| fixed | time (sine of radians) | -0.011 | -0.052 | 0.031 | ||
| fixed | time (cosine of radians) | 0.009 | -0.025 | 0.044 | ||
| fixed | temperaturre | 0.021 | -0.018 | 0.059 | ||
| fixed | Google Mobility | -0.09 | -0.222 | 0.045 | ||
| random | species within day & year (Intercept) | 5% | 5% | 6% | ||
| random | species within site (Intercept) | 4% | 4% | 5% | ||
| random | site (Intercept) | 8% | 8% | 8% | ||
| random | species (Intercept) | 16% | 14% | 17% | ||
| random | genus (Intercept) | 0% | 0% | 0% | ||
| random | Google Mobility (slope) | genus | 0% | 0% | 0% | ||
| random | weekday (Intercept) | 0% | 0% | 0% | ||
| random | country (Intercept) | 0% | -5% | 4% | ||
| random | Google Mobility (slope) | country | 0% | -5% | 4% | ||
| random | Residual | 66% | 56% | 77% | ||
| 0 2c | 3545 | fixed | (Intercept) | -0.056 | -0.141 | 0.028 |
| fixed | year | 0.014 | -0.024 | 0.053 | ||
| fixed | starting distance (ln) | 0.551 | 0.519 | 0.583 | ||
| fixed | flock size (ln) | -0.017 | -0.044 | 0.011 | ||
| fixed | body mass (ln) | -0.026 | -0.055 | 0.004 | ||
| fixed | time (sine of radians) | -0.033 | -0.074 | 0.009 | ||
| fixed | time (cosine of radians) | 0.018 | -0.018 | 0.053 | ||
| fixed | temperaturre | 0.031 | -0.008 | 0.068 | ||
| fixed | Google Mobility | -0.084 | -0.209 | 0.032 | ||
| random | site (Intercept) | 11% | 9% | 12% | ||
| random | country (Intercept) | 0% | -1% | 2% | ||
| random | Google Mobility (slope) | country | 0% | -1% | 2% | ||
| random | Residual | 89% | 84% | 93% | ||
| 0 3a | 3399 | fixed | (Intercept) | -0.024 | -0.178 | 0.131 |
| fixed | year | 0.027 | -0.017 | 0.071 | ||
| fixed | starting distance (ln) | 0.487 | 0.454 | 0.52 | ||
| fixed | flock size (ln) | 0.002 | -0.025 | 0.029 | ||
| fixed | body mass (ln) | 0.025 | -0.061 | 0.113 | ||
| fixed | time (sine of radians) | 0 | -0.041 | 0.045 | ||
| fixed | time (cosine of radians) | 0.007 | -0.028 | 0.043 | ||
| fixed | temperaturre | 0.022 | -0.018 | 0.06 | ||
| fixed | Google Mobility | -0.11 | -0.235 | 0.019 | ||
| random | species within day & year (Intercept) | 6% | 5% | 6% | ||
| random | species within site (Intercept) | 4% | 4% | 5% | ||
| random | site (Intercept) | 9% | 8% | 9% | ||
| random | species (Intercept) | 7% | 6% | 9% | ||
| random | genus (Intercept) | 3% | 3% | 4% | ||
| random | weekday (Intercept) | 0% | 0% | 0% | ||
| random | country (Intercept) | 0% | -5% | 4% | ||
| random | Google Mobility (slope) | country | 0% | -5% | 4% | ||
| random | Residual | 70% | 60% | 81% | ||
| 0 3b | 3399 | fixed | (Intercept) | -0.011 | -0.135 | 0.112 |
| fixed | year | 0.031 | -0.012 | 0.073 | ||
| fixed | starting distance (ln) | 0.487 | 0.454 | 0.518 | ||
| fixed | flock size (ln) | 0.003 | -0.024 | 0.03 | ||
| fixed | body mass (ln) | 0.036 | -0.048 | 0.122 | ||
| fixed | time (sine of radians) | -0.002 | -0.044 | 0.042 | ||
| fixed | time (cosine of radians) | 0.007 | -0.026 | 0.042 | ||
| fixed | temperaturre | 0.023 | -0.017 | 0.063 | ||
| fixed | Google Mobility | -0.128 | -0.258 | -0.003 | ||
| random | species within day & year (Intercept) | 6% | 5% | 6% | ||
| random | species within site (Intercept) | 4% | 4% | 5% | ||
| random | site (Intercept) | 9% | 9% | 9% | ||
| random | species (Intercept) | 8% | 7% | 9% | ||
| random | genus (Intercept) | 0% | -1% | 4% | ||
| random | Google Mobility (slope) | genus | 0% | -1% | 4% | ||
| random | weekday (Intercept) | 0% | 0% | 0% | ||
| random | country (Intercept) | 1% | -2% | 3% | ||
| random | Google Mobility (slope) | country | 1% | -2% | 3% | ||
| random | Residual | 69% | 58% | 79% | ||
| 0 3c | 3399 | fixed | (Intercept) | -0.068 | -0.155 | 0.02 |
| fixed | year | 0.023 | -0.019 | 0.063 | ||
| fixed | starting distance (ln) | 0.549 | 0.517 | 0.583 | ||
| fixed | flock size (ln) | -0.019 | -0.047 | 0.009 | ||
| fixed | body mass (ln) | -0.019 | -0.049 | 0.011 | ||
| fixed | time (sine of radians) | -0.024 | -0.066 | 0.018 | ||
| fixed | time (cosine of radians) | 0.015 | -0.021 | 0.053 | ||
| fixed | temperaturre | 0.03 | -0.01 | 0.069 | ||
| fixed | Google Mobility | -0.091 | -0.204 | 0.023 | ||
| random | site (Intercept) | 11% | 10% | 12% | ||
| random | country (Intercept) | 0% | 0% | 2% | ||
| random | Google Mobility (slope) | country | 0% | 0% | 2% | ||
| random | Residual | 88% | 83% | 91% |
Note that (1a) model is the one reported in the main text.
# TODO: modelAssCorrelations among predictors
d[, sin_rad := sin(rad)]
d[, cos_rad := cos(rad)]
dp <- d[, c("SD_ln", "flock_ln", "body_ln", "sin_rad", "cos_rad", "Temp", "Day")]
setnames(dp, old = c("SD_ln", "flock_ln", "body_ln", "sin_rad", "cos_rad", "Temp", "Day"), new = c("Starting distance\nln(m)", "Flock size\nln(m)", "Body mass\nln(m)", "Sine\nof radians", "Cosine\nof radians", "Temperature\n°C", "Day"))
#if (save_plot == TRUE) {
# png(here::here("Outputs/Fig_S1_rev.png"), width = 19, height = 19, units = "cm", bg = "transparent", res = 600)
# chart.Correlation(dp, histogram = TRUE, pch = 19, alpha = 0.5)
# mtext("Single observations", side = 3, line = 3)
# dev.off()
#}
chart.Correlation(dp, histogram = TRUE, pch = 19, alpha = 0.5)
mtext("Single observations", side = 3, line = 3)
Figure S2 | Pairwise correlations among fixed effects used in this study. On the diagonal: histograms and density lines (red) for each variable. Above diagonal: Pearson’s correlation with stars indicating significance . Below diagonal: the bivariate scatterplots, with each dot representing a single observation and a red line representing smoothed fit. Created with chart.Correlation function from R-package PerformanceAnalytics (Peterson & Carl 2020).
Species-site specific distributions
px = pp[N_during > 4 & N_before > 4]
dxx = d[paste(IDLocality, Species) %in% paste(px$IDLocality, px$Species)]
# table(dxx$IDLocality, dxx$Year)
#length(unique(px$IDLocality))
#length(unique(px$Species))
#sum(px$N_during) + sum(px$N_before)
dxx[, Country := factor(Country, levels = rev(c("Finland", "Poland", "Czechia", "Hungary", "Australia")))]
dxx[, sp_C_loc2 := paste(gsub("[_]", " ", Species), Country, IDLocality, sep = "\n")]
dxx[, genus := sub("_.*", "", Species)]
dxx[Covid == 0, period := "before COVID-19"]
dxx[Covid == 1, period := "during COVID-19"]
col3_ <- c("#357EBDFF", "#D43F3AFF", "#46B8DAFF", "#5CB85CFF", "#EEA236FF", "#9632B8FF", "#9632B8FF")[7:1]
col3__ <- col3_[3:7]
ggplot(dxx, aes(x = as.factor(Year), y = FID, col = Country, fill = period)) +
geom_boxplot(outlier.size = 0.5) +
facet_wrap(~sp_C_loc2) +
scale_y_continuous("Flight initiation distance [m]", trans = "log10") +
scale_x_discrete("Year", guide = guide_axis(angle = 45)) +
# scale_color_continuous() +
scale_colour_manual(values = col3__, guide = guide_legend(reverse = TRUE)) +
scale_fill_manual(values = c("white", "lightgrey")) +
theme_MB +
theme(
plot.title = element_text(size = 7),
strip.background = element_blank(),
strip.text.x = element_text(size = 5, color = "grey30", margin = margin(1, 1, 1, 1, "mm")),
# panel.spacing = unit(1, "mm"),
legend.position = "none", # c(1, 0.01),
legend.justification = c(1, 0),
legend.title = element_blank(),
# legend.spacing.y = unit(-0.78, "cm")
# legend.spacing.y = unit(0.02, "cm") use if LOESS smooth text as legend
legend.spacing.y = unit(-0.9, "cm"),
axis.text.x = element_text(colour = "grey30", size = 6),
axis.text.y = element_text(colour = "grey30", size = 6)
)
if(save_plot == TRUE){ggsave(here::here("Outputs/Fig_2_rev_v2.png"), width = 18, height = 16, units = "cm")}Figure 2 | Temporal variation in avian tolerance toward humans across species. Each heading denotes the scientific name of the species, country and unique site ID within each country. Boxplots outline colour highlights country (as in Fig. 1), fill colour indicates Period (white: before the COVID-19 shutdowns; grey: during the COVID-19 shutdowns). Boxplots depict median (horizontal line inside the box), the 25th and 75th percentiles (box) ± 1.5 times the interquartile range or the minimum/maximum value, whichever is smaller (bars), and the outliers (dots). Included are only species–site combinations with ≥5 observations per Period. Y-axis is on the log-scale. Note the lack of consistent shutdowns effects within and between species as well as within and between the countries.
Between/within-genus variation
# prepare data for fig 3 & s3
dxx <- d[paste(IDLocality, Species) %in% paste(pp$IDLocality, pp$Species)]
# length(dxx[, unique(paste(sp_loc, Covid))])
m <- lm(log(FID) ~ log(SD), dxx)
dxx[, resid_FID := resid(m)]
a <- dxx[, .(mean(resid_FID), sd(resid_FID), mean(FID), .N), by = .(Country, IDLocality, genus, Species, sp_loc, Covid)]
setnames(a, old = c("V1", "V2", "V3"), new = c("resid_FID_avg", "SD", "FID_avg"))
a[is.na(SD), SD := 0]
aw <- reshape(a, idvar = c("Country", "IDLocality", "genus", "Species", "sp_loc"), timevar = "Covid", direction = "wide")
aw[, Species := gsub("[_]", " ", Species)]
aw <- merge(aw, t, all.x = TRUE)
# table(aw$Family)
x <- aw[, .N, by = Species]
# x[order(Species)]
aw[, genus2 := genus]
aw[Species %in% x[N %in% c(1, 2), Species], genus2 := "other"]
# aw[genus2 == "Phoenicurus", unique(Species)]
ph[genus2 == "Motacilla" | uid %in% c("67a9ecfd-58ba-44a4-9986-243b6e610419"), uid := "cf522e02-35cc-44f5-841c-0e642987c2e4"]
ph[genus2 == "Sylvia", uid := "67a9ecfd-58ba-44a4-9986-243b6e610419"]
ph[, size := 0.2]
ph[genus2 %in% c("Anas", "Columba", "Dendrocopos", "Sturnus"), size := c(0.25, 0.25, 0.15, 0.1)]
ph[, FID_avg.0 := 1.5]
ph[, FID_avg.1 := 20]
ph[genus2 %in% c("Anas", "Columba"), FID_avg.0 := c(1.7, 1.7)]
ph[, resid_FID_avg.0 := -1.7]
ph[, resid_FID_avg.1 := 0.7]
ph[genus2 %in% c("Anas", "Columba"), resid_FID_avg.0 := c(-1.6, -1.6)]
ph[, genus2 := factor(genus2, levels = c("Anas", "Larus", "Columba", "Dendrocopos", "Picus", "Motacilla", "Erithacus", "Phoenicurus", "Turdus", "Sylvia", "Parus", "Sitta", "Pica", "Garrulus", "Corvus", "Sturnus", "Passer", "Fringilla", "other"))]
aw[, genus2 := factor(genus2, levels = c("Anas", "Larus", "Columba", "Dendrocopos", "Picus", "Motacilla", "Erithacus", "Phoenicurus", "Turdus", "Sylvia", "Parus", "Sitta", "Pica", "Garrulus", "Corvus", "Sturnus", "Passer", "Fringilla", "other"))]
aw[, genus2 := factor(genus2, levels = c("Anas", "Larus", "Columba", "Dendrocopos", "Picus", "Motacilla", "Erithacus", "Phoenicurus", "Turdus", "Sylvia", "Parus", "Sitta", "Pica", "Garrulus", "Corvus", "Sturnus", "Passer", "Fringilla", "other"))]
# Fig 3 and left panel of S3- plot from files
anas <- rasterGrob(change_col("#CCCCCC", readPNG(here::here("Data/Pics/Anas.png"))))
columba <- rasterGrob(change_col("#CCCCCC", readPNG(here::here("Data/Pics/Columba.png"))))
Dendrocopos <- rasterGrob(change_col("#CCCCCC", readPNG(here::here("Data/Pics/Dendrocopos.png"))))
Larus <- rasterGrob(change_col("#CCCCCC", readPNG(here::here("Data/Pics/Larus_flip.png"))))
Picus <- rasterGrob(change_col("#CCCCCC", readPNG(here::here("Data/Pics/Picus.png"))))
Motacilla <- rasterGrob(change_col("#CCCCCC", readPNG(here::here("Data/Pics/Motacilla.png"))))
Erithacus <- rasterGrob(change_col("#CCCCCC", readPNG(here::here("Data/Pics/Erithacus.png"))))
Phoenicurus <- rasterGrob(change_col("#CCCCCC", readPNG(here::here("Data/Pics/Phoenicurus.png"))))
Turdus <- rasterGrob(change_col("#CCCCCC", readPNG(here::here("Data/Pics/Turdus.png"))))
Sylvia <- rasterGrob(change_col("#CCCCCC", readPNG(here::here("Data/Pics/Sylvia.png"))))
Parus <- rasterGrob(change_col("#CCCCCC", readPNG(here::here("Data/Pics/Parus_flip.png"))))
Sitta <- rasterGrob(change_col("#CCCCCC", readPNG(here::here("Data/Pics/Sitta.png"))))
Pica <- rasterGrob(change_col("#CCCCCC", readPNG(here::here("Data/Pics/Pica.png"))))
Garrulus <- rasterGrob(change_col("#CCCCCC", readPNG(here::here("Data/Pics/Garrulus.png"))))
Corvus <- rasterGrob(change_col("#CCCCCC", readPNG(here::here("Data/Pics/Corvus.png"))))
Sturnus <- rasterGrob(change_col("#CCCCCC", readPNG(here::here("Data/Pics/Sturnus.png"))))
Passer <- rasterGrob(change_col("#CCCCCC", readPNG(here::here("Data/Pics/Passer_flip.png"))))
Fringilla <- rasterGrob(change_col("#CCCCCC", readPNG(here::here("Data/Pics/Fringilla.png"))))
other <- rasterGrob(change_col("#CCCCCC", readPNG(here::here("Data/Pics/other_flip.png"))))
ann_text <- data.frame(
FID_avg.0 = 8, FID_avg.1 = 10, lab = "Text",
genus2 = factor("Anas", levels = c("Anas", "Larus", "Columba", "Dendrocopos", "Picus", "Motacilla", "Erithacus", "Phoenicurus", "Turdus", "Sylvia", "Parus", "Sitta", "Pica", "Garrulus", "Corvus", "Sturnus", "Passer", "Fringilla", "other"))
)
ann_text2 <- data.frame(
FID_avg.0 = 6, FID_avg.1 = 3, lab = "Text",
genus2 = factor("Larus", levels = c("Anas", "Larus", "Columba", "Dendrocopos", "Picus", "Motacilla", "Erithacus", "Phoenicurus", "Turdus", "Sylvia", "Parus", "Sitta", "Pica", "Garrulus", "Corvus", "Sturnus", "Passer", "Fringilla", "other"))
)
aw2 <- data.frame(FID_avg.0 = c(11.25, 11.25), FID_avg.1 = c(3.5, 5.8), genus2 = factor("Larus", levels = c("Anas", "Larus", "Columba", "Dendrocopos", "Picus", "Motacilla", "Erithacus", "Phoenicurus", "Turdus", "Sylvia", "Parus", "Sitta", "Pica", "Garrulus", "Corvus", "Sturnus", "Passer", "Fringilla", "other")))
col3_ <- c("#357EBDFF", "#D43F3AFF", "#46B8DAFF", "#5CB85CFF", "#EEA236FF", "#9632B8FF", "#9632B8FF")[7:1]
col3__ <- col3_[3:7]
aw[, Country := factor(Country, levels = rev(c("Finland", "Poland", "Czechia", "Hungary", "Australia")))]g_gen <-
ggplot(aw, aes(x = FID_avg.0, y = FID_avg.1)) +
# geom_errorbar(aes(ymin = FID_avg.1-SD.1, ymax = FID_avg.1+SD.1, col = Country), width = 0) +
# geom_errorbar(aes(xmin = FID_avg.0-SD.0, xmax = FID_avg.0+SD.0, col = Country), width = 0) +
# geom_point(pch = 21, alpha = 0.7, aes(col = Country)) +
annotation_custom2(anas, data = ph[genus2 == "Anas"], xmin = 0.05, xmax = 0.5, ymax = 2.6) +
annotation_custom2(Larus, data = ph[genus2 == "Larus"], xmin = 0.05, xmax = 0.5, ymax = 2.6) +
annotation_custom2(columba, data = ph[genus2 == "Columba"], xmin = 0.05, xmax = 0.4, ymax = 2.7) +
annotation_custom2(Dendrocopos, data = ph[genus2 == "Dendrocopos"], xmin = 0.05, xmax = 0.25, ymax = 2.6) +
annotation_custom2(Picus, data = ph[genus2 == "Picus"], xmin = 0.05, xmax = 0.4, ymax = 2.7) +
annotation_custom2(Motacilla, data = ph[genus2 == "Motacilla"], xmin = 0.05, xmax = 0.5, ymax = 2.7) +
annotation_custom2(Erithacus, data = ph[genus2 == "Erithacus"], xmin = 0.05, xmax = 0.35, ymax = 2.7) +
annotation_custom2(Phoenicurus, data = ph[genus2 == "Phoenicurus"], xmin = 0.05, xmax = 0.35, ymax = 2.7) +
annotation_custom2(Turdus, data = ph[genus2 == "Turdus"], xmin = 0.05, xmax = 0.5, ymax = 2.7) +
annotation_custom2(Sylvia, data = ph[genus2 == "Sylvia"], xmin = 0.05, xmax = 0.5, ymax = 2.7) +
annotation_custom2(Parus, data = ph[genus2 == "Parus"], xmin = 0.05, xmax = 0.42, ymax = 2.7) +
annotation_custom2(Sitta, data = ph[genus2 == "Sitta"], xmin = 0.05, xmax = 0.5, ymax = 2.7) +
annotation_custom2(Pica, data = ph[genus2 == "Pica"], xmin = 0.05, xmax = 0.5, ymax = 2.5) +
annotation_custom2(Garrulus, data = ph[genus2 == "Garrulus"], xmin = 0.05, xmax = 0.6, ymax = 2.7) +
annotation_custom2(Corvus, data = ph[genus2 == "Corvus"], xmin = 0.05, xmax = 0.4, ymax = 2.55) +
annotation_custom2(Sturnus, data = ph[genus2 == "Sturnus"], xmin = 0.05, xmax = 0.24, ymax = 2.65) +
annotation_custom2(Passer, data = ph[genus2 == "Passer"], xmin = 0.05, xmax = 0.36, ymax = 2.65) +
annotation_custom2(Fringilla, data = ph[genus2 == "Fringilla"], xmin = 0.05, xmax = 0.5, ymax = 2.7) +
annotation_custom2(other, data = ph[genus2 == "other"], xmin = 0.05, xmax = 0.45, ymax = 2.4) +
geom_point(pch = 21, alpha = 0.7, aes(fill = Country), col = "grey20", alpha = 0.8) + # col = "white") +
# ggtitle ("Sim based")+
geom_abline(intercept = 0, slope = 1, lty = 3, col = "grey80") +
geom_line(data = aw2, col = "grey80", lwd = 0.25) +
geom_text(data = ann_text, label = "No difference", col = "grey80", angle = 45, size = 2) +
geom_text(data = ann_text2, label = "Species mean / site", col = "grey60", size = 2, ) +
facet_wrap(~genus2) +
# geom_phylopic(data = o, aes(image = uid), color = "grey80", size = o$size) + # ,
# scale_fill_viridis(discrete = TRUE, guide = guide_legend(reverse = FALSE)) +
scale_fill_manual(values = col3__, guide = guide_legend(reverse = TRUE)) +
scale_x_continuous("Before COVID-19 shutdown - flight initiation distance [m]", expand = c(0, 0), trans = "log10") +
scale_y_continuous("During COVID-19 shutdown - flight initiation distance [m]", expand = c(0, 0), trans = "log10") +
# labs(title = "Species means per sampling location")+
theme_MB +
theme(
plot.title = element_text(size = 7),
strip.background = element_blank(),
# panel.spacing = unit(1, "mm"),
legend.position = c(1, 0.025),
legend.justification = c(1, -0.05)
)
gg_gen <- ggplotGrob(g_gen) # gg$layout$name
ggx_gen <- gtable_filter_remove(gg_gen,
name = paste0("axis-b-", c(2, 4), "-4"),
trim = FALSE
)
if (save_plot == TRUE) {
ggsave(here::here("Outputs/Fig_3_width-122mm_col_grey_rev.png"), ggx_gen, width = 4.8, height = 4.5, dpi = 600) # 12.2cm # with label inside
}
grid.draw(ggx_gen)
Figure 3 | Avian tolerance towards humans before and during the COVID-19 shutdowns according to genera. Dots represent means or single escape distance observations of species at specific sites (e.g. park or cemetery) with data for both periods (before and during the shutdowns) and not corrected for other factors such as starting distance of the observer (plot corrected for starting distance gives similar patterns: Fig. S3). Dot colour highlights the country. Dotted lines indicate no difference; dots above the lines indicate lower tolerance towards humans (i.e. longer escape distances), dots below the lines indicate higher tolerance during than before the COVID-19 shutdowns. Panels are ordered according to evolutionary history of birds with top left panels representing the oldest genera, and bottom right, the youngest. Panel titled ‘other’ contains genera with only one or two data points. The axes are on the log-scale. For a species-specific figure, see Fig. S4. Silhouette of Garrulus glandarius, Motacilla alba, Picus viridis, Phoenicurus ochruros, Sylvia borin were drawn by Martin Bulla, Erithacus rubecula drawn by Rebecca Groom, and Fringilla coelebs and Sturnus vulgaris by Maxime Dahirel and all are available at PhyloPic under Creative Commons Attribution 3.0 Unported licence. The remaining silhouettes are available at PhyloPic under the Public Domain Dedication 1.0 license.
# Fig S3 right panel
g_2 <-
ggplot(aw, aes(x = resid_FID_avg.0, y = resid_FID_avg.1)) +
geom_point(pch = 21, alpha = 0.7, aes(fill = Country), col = "grey20", alpha = 0.8) + # col = "white") +
geom_abline(intercept = 0, slope = 1, lty = 3, col = "grey80") +
facet_wrap(~genus2) +
# geom_phylopic(data = o, aes(image = uid), color = "grey80", size = o$size) +
# scale_fill_viridis(discrete = TRUE, guide = guide_legend(reverse = TRUE)) +
scale_fill_manual(values = col3__, guide = guide_legend(reverse = TRUE)) +
scale_x_continuous("Before COVID-19 shutdown - residual escape distance", expand = c(0, 0)) +
scale_y_continuous("During COVID-19 shutdown - residual escape distance", expand = c(0, 0)) +
# labs(title = "Species means per sampling location")+
theme_MB +
theme(
plot.title = element_text(size = 7),
strip.background = element_blank(),
legend.position = "none",
# legend.position = c(1, 0),
legend.justification = c(1, 0)
)
g_g2 <- ggplotGrob(g_2) # gg$layout$name
g_gx2 <- gtable_filter_remove(g_g2,
name = paste0("axis-b-", c(2, 4), "-4"),
trim = FALSE
)
# Fig S5 combine
grid.draw(cbind(ggx_gen, g_gx2, size = "last"))
if (save_plot == TRUE) {
ggsave(here::here("Outputs/Fig_S3_rev_v3.png"), cbind(ggx_gen, g_gx2, size = "last"), width = 4.8 * 2, height = 4.5, dpi = 600)
}Figure S3 | Comparison of genus-specific flight initiation distance (left) and residual flight initiation distance (right) before and during the COVID-19 shutdown. Left panel is a copy of a main text Fig. 3 (see there for details). The escape distance represents the raw data that can be confounded by the observers starting distance (for our data the rPearson = 0.58. Right panel dots depict residual escape distances from a model with flight initiation distance (ln-transformed) as a response and starting distance (ln-transformed) as a predictor, i.e. the dots represent before and during shutdowns values that are controlled for starting distance. Note that such control for starting distance (right) has little influence on the depicted relationships. Indeed, the Pearson’s correlation coefficient for escape distance (ln-scale) and residual escape distance was 0.8 for single values and 0.74 for the species means per sampling location.
Difference in number of observations per species and site before and during shutdowns:
ggplot(aw, aes(x = N.0 - N.1)) +
geom_histogram() + xlab('Before minus during shutdowns\n[# observations')
# nrow(aw[abs(N.0 - N.1) > 2])
# nrow(aw[!abs(N.0 - N.1) > 2])Between/within-species variation
aw[, sp2 := gsub(" ", "\n", Species)]
ann_text <- data.frame(
FID_avg.0 = 8, FID_avg.1 = 10, lab = "Text",
Species = factor("Aegithalos caudatus", levels = levels(as.factor(aw$Species)))
)
ann_text$sp2 = gsub(" ", "\n", ann_text$Species)
g3 <-
ggplot(aw, aes(x = FID_avg.0, y = FID_avg.1)) +
# geom_errorbar(aes(ymin = FID_avg.1-SD.1, ymax = FID_avg.1+SD.1, col = Country), width = 0) +
# geom_errorbar(aes(xmin = FID_avg.0-SD.0, xmax = FID_avg.0+SD.0, col = Country), width = 0) +
# geom_point(pch = 21, alpha = 0.7, aes(col = Country)) +
geom_point(pch = 21, alpha = 0.7, aes(fill = Country), col = "grey20", alpha = 0.8) + # col = "white") +
# ggtitle ("Sim based")+
geom_abline(intercept = 0, slope = 1, lty = 3, col = "grey80") +
geom_text(data = ann_text, label = "No difference", col = "grey80", angle = 45, size = 2) +
facet_wrap(~sp2) +
# geom_phylopic(data = o, aes(image = uid), color = "grey80", size = o$size) + # ,
# scale_fill_viridis(discrete=TRUE,guide = guide_legend(reverse = FALSE)) +
scale_fill_manual(values = col3__, guide = guide_legend(reverse = TRUE)) +
scale_x_continuous("Before COVID-19 shutdown - flight initiation distance [m]", expand = c(0, 0), trans = "log10") +
scale_y_continuous("During COVID-19 shutdown - flight initiation distance [m]", expand = c(0, 0), trans = "log10") +
labs(title = "Species means per sampling location") +
theme_MB +
theme(
plot.title = element_text(size = 7),
strip.background = element_blank(),
# panel.spacing = unit(1, "mm"),
legend.position = c(0.96, 0.0),
legend.justification = c(1, 0)
)
gg3 <- ggplotGrob(g3) # gg2$layout$name
ggx3 <- gtable_filter_remove(gg3,
name = c(paste0("axis-b-", c(2, 4), "-7"), "axis-b-6-6"),
trim = FALSE
)
grid.draw(ggx3)
if (save_plot == TRUE) {
ggsave(here::here("Outputs/Fig_S4_species_rev_v4.png"), ggx3, width = 13.5, height = 17.5, unit = "cm", dpi = 600) # 11.43cm
}Figure S4 | Avian tolerance towards humans before and during the COVID-19 shutdowns according to species. Dots represent means or single escape distance observations of species at specific sites (e.g. specific park or cemetery) with data for both periods (i.e. before and during the shutdowns) and not corrected for other factors such as starting distance of the observer. Dot colour highlights the country. Dotted lines indicate no difference, dots above the lines indicate lower tolerance towards humans (i.e. longer escape distances), dots below the lines indicate lower tolerance before than during the COVID-19 shutdowns. Panels are ordered alphabetically. The axes are on the log-scale.
Exploration of Google Mobility
g_ <- fread(here::here("Data/google_mobility.txt")) # fwrite(d, here::here('Data/data.txt'), sep ='\t')
g_[, Year := as.integer(substring(date, nchar(date) - 3, nchar(date)))]
g_[nchar(date) == 9, date := paste0("0", date)]
g_[, date_ := as.Date(date, format = "%d.%m.%Y")]
g_[, Day := yday(date_)]
setnames(g_, old = "country_region", new = "Country")
g_[, weekday := weekdays(date_)]
g0 = ggplot(g_, aes(x = parks_percent_change_from_baseline, fill = factor(Year))) +
geom_histogram(position = "dodge") +
# scale_y_continuous(trans = 'log')+
scale_fill_manual(values = c("orange", "skyblue", "black"), guide = 'none') +
geom_vline(xintercept = 0, lty = 3, col = "red") +
labs(subtitle = "Distribution") +
xlab( 'Google Mobility\n[% change in human presence]') +
ylab( 'Count') +
facet_wrap(~Country, nrow = 5)
g1 = ggplot(g_, aes(x = Day, y = parks_percent_change_from_baseline, col = factor(Year))) +
geom_line() +
facet_wrap(~Country, nrow = 5) +
labs(subtitle = "Raw data", xlab = 'Day\n ', y = 'Google Mobility\n[% change in human presence]') +
# scale_y_continuous(trans = 'log')+
coord_cartesian(ylim = c(-100, 300))+
scale_color_manual(values = c("orange", "skyblue", "black"), guide = 'none')
g2 = ggplot(g_, aes(x = Day, y = parks_percent_change_from_baseline, col = factor(Year))) +
stat_smooth() +
facet_wrap(~Country, nrow = 5) +
labs(subtitle = "Loess", xlab = 'Day\n ') +
# scale_y_continuous(trans = 'log')+
coord_cartesian(ylim = c(-100, 300))+
scale_color_manual(values = c("orange", "skyblue", "black"), name = 'Year')+
theme(
axis.text.y = element_blank(),
axis.title.y = element_blank()
)
ggarrange(
g0, g1, g2,
ncol = 3, widths = c(0.9, 1, 1.1)
)
ggsave(here::here("Outputs/Fig_4_rev.png"), width = 8*2.54, height = 5*2.54, unit = "cm", dpi = 600)Figure 4 | Changes in human presence in parks within and between years and countries. Left plots represent distribution (histograms) of human presence (Google Mobility), middle plots the raw data, right plots locally estimated scatterplot smoothing. Note that Google Mobility data were not freely available for years before the COVID-19 pandemic (i.e. before 2020) but 2022 was a year without shutdowns in the studied countries. For weekday-specific pattern see Fig. S5 below.
g[, weekday := factor(weekday, levels = (c("Monday", "Tuesday", "Wednesday", "Thursday", "Friday", "Saturday", "Sunday")))]
ggplot(g, aes(x = Day, y = parks_percent_change_from_baseline, col = factor(Year))) +
geom_line() +
facet_grid(rows = vars(Country), cols = vars(weekday)) +
# scale_y_continuous(trans = 'log')+
scale_color_manual(values = c("orange", "skyblue", "black"))
ggsave(here::here("Outputs/Fig_S5_rev.png"), width = 8 * 2.54, height = 6 * 2.54, unit = "cm", dpi = 600)Figure S5 | Changes in human presence (Google Mobility) in parks across weekdays and years. Depicted are raw data connected by lines. Note that Google Mobility data were not freely available for years before the COVID-19 pandemic (i.e. before 2020) but 2022 was a year without shutdowns in the studied countries.
Google Mobility ~ stringency
# Predictions
l = list()
sc = s[Country == "Czechia"]
cz <- lmer(parks_percent_change_from_baseline ~
StringencyIndex +
(scale(StringencyIndex)|weekday),
# (1 | Year) + (1 | weekday) + (1|genus) + (1 | Species) + (1 | sp_day_year) + (1 | IDLocality),
data = sc, REML = FALSE
)
bsim <- sim(cz, n.sim = nsim)
v <- apply(bsim@fixef, 2, quantile, prob = c(0.5))
newD <- data.frame(StringencyIndex = seq(min(sc$StringencyIndex), max(sc$StringencyIndex), length.out = 100)) # values to predict for
X <- model.matrix(~StringencyIndex, data = newD) # exactly the model which was used has to be specified here
newD$pred <- (X %*% v)
predmatrix <- matrix(nrow = nrow(newD), ncol = nsim)
for (j in 1:nsim) predmatrix[, j] <- (X %*% bsim@fixef[j, ])
newD$lwr <- apply(predmatrix, 1, quantile, prob = 0.025)
newD$upr <- apply(predmatrix, 1, quantile, prob = 0.975)
newD$pred <- apply(predmatrix, 1, quantile, prob = 0.5)
newD$Country = 'Czechia'
l[[1]] = newD
s[, year_weekday :=paste(Year, weekday)]
sf = s[Country == "Finland"]
fi <- lmer(parks_percent_change_from_baseline ~
Year+
StringencyIndex +
(scale(StringencyIndex)|year_weekday),
# (1 | Year) + (1 | weekday) + (1|genus) + (1 | Species) + (1 | sp_day_year) + (1 | IDLocality),
data = sf, REML = FALSE
)
bsim <- sim(fi, n.sim = nsim)
v <- apply(bsim@fixef, 2, quantile, prob = c(0.5))
newD <- data.frame(Year = mean(sf$Year), StringencyIndex = seq(min(sf$StringencyIndex), max(sf$StringencyIndex), length.out = 100)) # values to predict for
X <- model.matrix(~Year + StringencyIndex, data = newD) # exactly the model which was used has to be specified here
newD$pred <- (X %*% v)
predmatrix <- matrix(nrow = nrow(newD), ncol = nsim)
for (j in 1:nsim) predmatrix[, j] <- (X %*% bsim@fixef[j, ])
newD$lwr <- apply(predmatrix, 1, quantile, prob = 0.025)
newD$upr <- apply(predmatrix, 1, quantile, prob = 0.975)
newD$pred <- apply(predmatrix, 1, quantile, prob = 0.5)
newD$Country = 'Finland'
newD$Year = NULL
l[[2]] = newD
sh <- s[Country == "Hungary"]
hu <- lmer(parks_percent_change_from_baseline ~
Year +
StringencyIndex +
(scale(StringencyIndex) | year_weekday),
# (1 | Year) + (1 | weekday) + (1|genus) + (1 | Species) + (1 | sp_day_year) + (1 | IDLocality),
data = sh, REML = FALSE
)
bsim <- sim(hu, n.sim = nsim)
v <- apply(bsim@fixef, 2, quantile, prob = c(0.5))
newD <- data.frame(Year = mean(sh$Year), StringencyIndex = seq(min(sh$StringencyIndex), max(sh$StringencyIndex), length.out = 100)) # values to predict for
X <- model.matrix(~ Year + StringencyIndex, data = newD) # exactly the model which was used has to be specified here
newD$pred <- (X %*% v)
predmatrix <- matrix(nrow = nrow(newD), ncol = nsim)
for (j in 1:nsim) predmatrix[, j] <- (X %*% bsim@fixef[j, ])
newD$lwr <- apply(predmatrix, 1, quantile, prob = 0.025)
newD$upr <- apply(predmatrix, 1, quantile, prob = 0.975)
newD$pred <- apply(predmatrix, 1, quantile, prob = 0.5)
newD$Country <- "Hungary"
newD$Year <- NULL
l[[3]] <- newD
sp <- s[Country == "Poland"]
pl <- lmer(parks_percent_change_from_baseline ~
Year +
StringencyIndex +
(scale(StringencyIndex) | year_weekday),
# (1 | Year) + (1 | weekday) + (1|genus) + (1 | Species) + (1 | sp_day_year) + (1 | IDLocality),
data = sp, REML = FALSE
)
bsim <- sim(pl, n.sim = nsim)
v <- apply(bsim@fixef, 2, quantile, prob = c(0.5))
newD <- data.frame(Year = mean(sp$Year), StringencyIndex = seq(min(sp$StringencyIndex), max(sp$StringencyIndex), length.out = 100)) # values to predict for
X <- model.matrix(~ Year + StringencyIndex, data = newD) # exactly the model which was used has to be specified here
newD$pred <- (X %*% v)
predmatrix <- matrix(nrow = nrow(newD), ncol = nsim)
for (j in 1:nsim) predmatrix[, j] <- (X %*% bsim@fixef[j, ])
newD$lwr <- apply(predmatrix, 1, quantile, prob = 0.025)
newD$upr <- apply(predmatrix, 1, quantile, prob = 0.975)
newD$pred <- apply(predmatrix, 1, quantile, prob = 0.5)
newD$Country <- "Poland"
newD$Year <- NULL
l[[4]] <- newD
sa <- s[Country == "Australia"]
au <- lmer(parks_percent_change_from_baseline ~
Year +
StringencyIndex +
(scale(StringencyIndex) | year_weekday),
# (1 | Year) + (1 | weekday) + (1|genus) + (1 | Species) + (1 | sp_day_year) + (1 | IDLocality),
data = sa, REML = FALSE
)
bsim <- sim(au, n.sim = nsim)
v <- apply(bsim@fixef, 2, quantile, prob = c(0.5))
newD <- data.frame(Year = mean(sa$Year), StringencyIndex = seq(min(sa$StringencyIndex), max(sa$StringencyIndex), length.out = 100)) # values to predict for
X <- model.matrix(~ Year + StringencyIndex, data = newD) # exactly the model which was used has to be specified here
newD$pred <- (X %*% v)
predmatrix <- matrix(nrow = nrow(newD), ncol = nsim)
for (j in 1:nsim) predmatrix[, j] <- (X %*% bsim@fixef[j, ])
newD$lwr <- apply(predmatrix, 1, quantile, prob = 0.025)
newD$upr <- apply(predmatrix, 1, quantile, prob = 0.975)
newD$pred <- apply(predmatrix, 1, quantile, prob = 0.5)
newD$Country <- "Australia"
newD$Year <- NULL
l[[5]] <- newD
# Figure G_S
g_s = data.table(do.call(rbind,l))
g_s[, Country := factor(Country, levels = rev(c("Finland", "Poland", "Czechia", "Hungary", "Australia")))]
col3_ = c("#357EBDFF", "#D43F3AFF", "#46B8DAFF", "#5CB85CFF", "#EEA236FF", "#9632B8FF", "#9632B8FF")[7:1]
col3__ = col3_[3:7]
p =
ggplot(g_s, aes(x = StringencyIndex, y = pred, col = Country)) +
geom_ribbon(aes(ymin = lwr, ymax = upr, fill = Country, color = NULL), alpha = .15) +
geom_jitter(aes(y = parks_percent_change_from_baseline, fill = Country), data = s, pch = 21, col = 'grey20', width = 0.7, height = 3, alpha = 0.5) +
geom_line(lwd = 1) +
labs(subtitle = "Mixed model per country predicitons", y = "Google Mobiligy\n[% change from baseline]", x = "Stringency Index") +
# scale_color_locuszoom()+
# scale_fill_locuszoom(guide = "none")
scale_x_continuous(breaks = round(seq(25, 75, by = 25), 1)) +
scale_y_continuous(breaks = round(seq(-100, 200, by = 50), 1)) +
#scale_y_continuous(breaks = round(seq(-100, 175, by = 25), 1)) +
scale_colour_manual(values = col3__, guide = guide_legend(reverse = TRUE, override.aes = list(size = 0)),
labels = paste("<span style='color:",
col3__,
"'>",
levels(g_s$Country),
"</span>")
) +
scale_fill_manual(values = col3__, guide = "none") +
theme_bw() +
theme(
legend.text = element_markdown(size = 6),
#legend.position = "right",
legend.title = element_blank(),
# legend.spacing.y = unit(0.1, 'cm'),
legend.key.height = unit(0.5, "line"),
legend.key.size = unit(0, "line"),
legend.margin = margin(0, 0, 0, 0),
legend.box.margin=margin(-10,1,-10,-10),
# legend.position=c(0.5,1.6),
plot.title = element_text(color = "grey", size = 7),
plot.subtitle = element_text(color = "grey60", size = 6),
plot.margin = margin(b = 0.5, l = 0.5, t = 0.5, r = 0.5, unit = "pt"),
panel.grid = element_blank(),
panel.border = element_blank(),
panel.background = element_blank(),
axis.line = element_line(colour = ax_lines, size = 0.25),
axis.ticks = element_line(colour = ax_lines, size = 0.25),
# axis.text.x = element_text()
axis.ticks.length = unit(1, "pt"),
axis.text = element_text(, size = 6),
axis.title = element_text(size = 7)
)
p 
if(save_plot==TRUE){
ggsave(here::here("Outputs/Fig_5_rev.png"), p + theme(plot.subtitle = element_blank()), width = 7, height = 6, unit = "cm", dpi = 600)
}Figure 5 | Association between human presence in parks (Google Mobility) and stringency of antipandemic governmental restrictions (stringency index). Lines with shaded areas represent predicted relationships from country-specific mixed effect models controlled for the year and non-independence of data points by including weekday within the year as random intercept and stringency index as a random slope (Table S3). Dots represent raw data, jittered to increase visibility. Colours indicate country. Note the generally negative but weak association between human presence and stringency index.
Table S3 | Google Mobility in relation to stringency index
ll = list()
s[, year_weekday := paste(Year, weekday)]
sf = s[Country == "Finland"]
fi <- lmer(scale(parks_percent_change_from_baseline) ~
scale(Year) +
scale(StringencyIndex) +
(scale(StringencyIndex) | year_weekday),
data = sf, REML = FALSE
)
ll[[1]] = m_out(name = "Table S3 - FI", dep = "Google Mobility", model = fi, nsim = 5000)
sp <- s[Country == "Poland"]
pl <- lmer(scale(parks_percent_change_from_baseline) ~
scale(Year) +
scale(StringencyIndex) +
(scale(StringencyIndex) | year_weekday),
data = sp, REML = FALSE
)
ll[[2]] = m_out(name = "Table S3 - PL", dep = "Google Mobility", model = pl, nsim = 5000)
sc = s[Country == "Czechia"]
cz <- lmer(scale(parks_percent_change_from_baseline) ~
scale(StringencyIndex) +
(scale(StringencyIndex) | weekday),
data = sc, REML = FALSE
)
ll[[3]] = m_out(name = "Table S3 - CZ", dep = "Google Mobility", model = cz, nsim = 5000)
sh <- s[Country == "Hungary"]
hu <- lmer(scale(parks_percent_change_from_baseline) ~
scale(Year) +
scale(StringencyIndex) +
(scale(StringencyIndex) | year_weekday),
data = sh, REML = FALSE
)
ll[[4]] = m_out(name = "Table S3 - HU", dep = "Google Mobility", model = hu, nsim = 5000)
sa <- s[Country == "Australia"]
au <- lmer(scale(parks_percent_change_from_baseline) ~
scale(Year) +
scale(StringencyIndex) +
(scale(StringencyIndex) | year_weekday),
# (1 | Year) + (1 | weekday) + (1|genus) + (1 | Species) + (1 | sp_day_year) + (1 | IDLocality),
data = sa, REML = FALSE
)
ll[[5]] = m_out(name = "Table S3 - AU", dep = "Google Mobility", model = au, nsim = 5000)
out_g_s = data.table(do.call(rbind, ll))
out_g_s[is.na(out_g_s)] <- ""
out_g_s$R2_mar = out_g_s$R2_con = NULL
out_g_s[, effect := gsub("scale\\(Year\\)", "year", effect)]
out_g_s[, effect := gsub("scale\\(StringencyIndex\\)", "stringency index", effect)]
out_g_s[, effect := gsub("year_weekday", "weekday within year", effect)]
out_g_s[type == "random" & grepl("stringency index", effect, fixed = TRUE), effect := paste("stringency index (slope) |", gsub(" stringency index", "", effect))]
fwrite(file = here::here("Outputs/Table_S3_rev.csv"), out_g_s)
out_g_s$error_structure = out_g_s$response = NULL
out_g_s[model!="", model:=c('Finland',
'Poland',
'Czechia',
'Hungary',
'Australia')]
setnames(out_g_s, old = c("estimate_r", "lwr_r", "upr_r"), new = c("estimate", "lower", "upper"))
out_g_s %>%
kbl() %>%
kable_paper("hover", full_width = F)| model | N | type | effect | estimate | lower | upper |
|---|---|---|---|---|---|---|
| Finland | 322 | fixed | (Intercept) | -0.27 | -1.02 | 0.529 |
| fixed | year | 0.936 | 0.441 | 1.44 | ||
| fixed | stringency index | -0.261 | -1.023 | 0.492 | ||
| random | weekday within year (Intercept) | 22% | 47% | 54% | ||
| random | stringency index (slope) | weekday within year | 22% | 47% | 54% | ||
| random | Residual | 55% | -7% | 5% | ||
| Poland | 329 | fixed | (Intercept) | 0.023 | -0.627 | 0.626 |
| fixed | year | 0.76 | 0.015 | 1.531 | ||
| fixed | stringency index | 0.016 | -0.359 | 0.4 | ||
| random | weekday within year (Intercept) | 47% | 44% | 49% | ||
| random | stringency index (slope) | weekday within year | 47% | 44% | 49% | ||
| random | Residual | 6% | 2% | 13% | ||
| Czechia | 1054 | fixed | (Intercept) | 0.009 | -0.521 | 0.548 |
| fixed | stringency index | -0.229 | -0.642 | 0.191 | ||
| random | weekday (Intercept) | 23% | -182% | 40% | ||
| random | stringency index (slope) | weekday | 23% | -182% | 40% | ||
| random | Residual | 54% | 21% | 464% | ||
| Hungary | 874 | fixed | (Intercept) | 0.249 | -0.255 | 0.766 |
| fixed | year | 0.002 | -0.442 | 0.465 | ||
| fixed | stringency index | -1.269 | -1.969 | -0.574 | ||
| random | weekday within year (Intercept) | 48% | 21% | 49% | ||
| random | stringency index (slope) | weekday within year | 48% | 21% | 49% | ||
| random | Residual | 4% | 2% | 59% | ||
| Australia | 1065 | fixed | (Intercept) | 0.519 | -0.011 | 1.053 |
| fixed | year | 0.443 | 0.241 | 0.645 | ||
| fixed | stringency index | -0.537 | -0.848 | -0.246 | ||
| random | weekday within year (Intercept) | -31% | 44% | 65% | ||
| random | stringency index (slope) | weekday within year | -31% | 44% | 65% | ||
| random | Residual | 162% | -30% | 13% |
Effect sizes for stringency or Google Mobility
# predictions for fig and table for stringency
# full
mss <- lmer(scale(log(FID)) ~
scale(Year) +
scale(log(SD)) +
scale(log(FlockSize)) +
scale(log(BodyMass)) +
scale(sin(rad)) + scale(cos(rad)) +
# scale(Day)+
scale(Temp) +
scale(StringencyIndex) +
(1|weekday) + (1 | genus) + (1 | Species) + (1 | sp_day_year) + (scale(StringencyIndex) | Country) + (1 | IDLocality) + (1 | sp_loc),
data = s, REML = FALSE,
control = lmerControl(
optimizer = "optimx", optCtrl = list(method = "nlminb")
)
)
est_mss <- est_out(mss, "ALL: (1|weekday) + (1|genus) + (1|Species) + (1|sp_day_year) + (scale(StringencyIndex)|Country) + (1|IDLocality) +(1|sp_loc)")
est_mss[, control_for_starting_distance := "yes"]
msx <- lmer(scale(log(FID)) ~
scale(Year) +
#scale(log(SD)) +
scale(log(FlockSize)) +
scale(log(BodyMass)) +
scale(sin(rad)) + scale(cos(rad)) +
# scale(Day)+
scale(Temp) +
scale(StringencyIndex) +
(1|weekday) + (1 | genus) + (1 | Species) + (1 | sp_day_year) + (scale(StringencyIndex) | Country) + (1| IDLocality) + (1 | sp_loc),
data = s, REML = FALSE,
control = lmerControl(
optimizer = "optimx", optCtrl = list(method = "nlminb")
)
)
est_msx <- est_out(msx, "ALL: (1|weekday) + (1|genus) + (1|Species) + (1|sp_day_year) + (scale(StringencyIndex)|Country) + (1|IDLocality) +(1|sp_loc)")
est_msx[, control_for_starting_distance := "no"]
# CZ - singular fits only due to genera estimated as zero (removing it changes no results)
css <- lmer(scale(log(FID)) ~
#scale(Year) +
scale(log(SD)) +
scale(log(FlockSize)) +
scale(log(BodyMass)) +
scale(sin(rad)) + scale(cos(rad)) +
# scale(Day)+
scale(Temp) +
scale(StringencyIndex) +
(1|weekday) + (1 | genus) + (1| Species) + (1 | sp_day_year) + (1| IDLocality) + (1|sp_loc),
#(1 | Year) + (1 | weekday) + (1|genus) + (1 | Species) + (1 | sp_day_year) + (1 | IDLocality),
data = s[Country == "Czechia"], REML = FALSE
)
est_css <- est_out(css, "Czechia: (1|weekday) + (1|genus)+(1|Species)+(1|sp_day_year)+(1|IDLocality)+(1|sp_loc)")
est_css[, control_for_starting_distance := "yes"]
csx <- lmer(scale(log(FID)) ~
#scale(Year) +
scale(log(SD)) +
scale(log(FlockSize)) +
scale(log(BodyMass)) +
scale(sin(rad)) + scale(cos(rad)) +
# scale(Day)+
scale(Temp) +
scale(StringencyIndex) +
(1|weekday) + (1 | genus) + (1| Species) + (1 | sp_day_year) + (1| IDLocality) + (1|sp_loc),
#(1 | Year) + (1 | weekday) + (1|genus) + (1 | Species) + (1 | sp_day_year) + (1 | IDLocality),
data = s[Country == "Czechia"], REML = FALSE
)
est_csx <- est_out(csx, "Czechia: (1|weekday) + (1|genus)+(1|Species)+(1|sp_day_year)+(1|IDLocality)+(1|sp_loc)")
est_csx[, control_for_starting_distance := "no"]
# FI
fss <- lmer(scale(log(FID)) ~
scale(Year) +
scale(log(SD)) +
scale(log(FlockSize)) +
scale(log(BodyMass)) +
scale(sin(rad)) + scale(cos(rad)) +
# scale(Day)+
scale(Temp) +
scale(StringencyIndex) +
(1|weekday) + (1 | genus) + (1| Species) + (1 | sp_day_year) + (1| IDLocality) + (1|sp_loc),
data = s[Country == "Finland"], REML = FALSE
)
est_fss <- est_out(fss, "Finland: (1|weekday) + (1|genus)+(1|Species)+(1|sp_day_year)+(scale(StringencyIndex)|IDLocality)+(1|sp_loc)")
est_fss[, control_for_starting_distance := "yes"]
fsx <- lmer(scale(log(FID)) ~
scale(Year) +
# scale(log(SD)) +
scale(log(FlockSize)) +
scale(log(BodyMass)) +
scale(sin(rad)) + scale(cos(rad)) +
# scale(Day)+
scale(Temp) +
scale(StringencyIndex) +
(1|weekday) + (1 | genus) + (1| Species) + (1 | sp_day_year) + (1| IDLocality) + (1|sp_loc),
data = s[Country == "Finland"], REML = FALSE
)
est_fsx <- est_out(fsx, "Finland: (1|weekday) + (1|genus)+(1|Species)+(1|sp_day_year)+(scale(StringencyIndex)|IDLocality)+(1|sp_loc)")
est_fsx[, control_for_starting_distance := "no"]
# HU - singular fits only due to sp_loc estimated as zero (removing it changes no results)
hss <- lmer(scale(log(FID)) ~
scale(Year) +
scale(log(SD)) +
scale(log(FlockSize)) +
scale(log(BodyMass)) +
scale(sin(rad)) + scale(cos(rad)) +
# scale(Day)+
scale(Temp) +
scale(StringencyIndex) +
(1|weekday) + (1 | genus) + (1 | Species) + (1 | sp_day_year) + (1 | IDLocality) + (1 | sp_loc),
# (1 | Year) + (1 | weekday) + (1|genus) + (1 | Species) + (1 | sp_day_year) + (1 | IDLocality),
data = s[Country == "Hungary"], REML = FALSE
)
est_hss <- est_out(hss, "Hungary: (1|weekday) + (1|genus)+(1|Species)+(1|sp_day_year)+(1|IDLocality)+(1|sp_loc)")
est_hss[, control_for_starting_distance := "yes"]
hsx <- lmer(scale(log(FID)) ~
scale(Year) +
# scale(log(SD)) +
scale(log(FlockSize)) +
scale(log(BodyMass)) +
scale(sin(rad)) + scale(cos(rad)) +
# scale(Day)+
scale(Temp) +
scale(StringencyIndex) +
(1|weekday) + (1 | genus) + (1 | Species) + (1 | sp_day_year) + (1 | IDLocality) + (1 | sp_loc),
# (1 | Year) + (1 | weekday) + (1 | genus) + (1 | Species) + (1 | sp_day_year) + (1 | IDLocality),
data = s[Country == "Hungary"], REML = FALSE
)
est_hsx <- est_out(hsx, "Hungary: (1|weekday) + (1|genus)+(1|Species)+(1|sp_day_year)+(1|IDLocality)+(1|sp_loc)")
est_hsx[, control_for_starting_distance := "no"]
# AU - singular fits only due to Year and random slope estimated as zero (removing those changes no results)
ass <- lmer(scale(log(FID)) ~
scale(Year) +
scale(log(SD)) +
scale(log(FlockSize)) +
scale(log(BodyMass)) +
scale(sin(rad)) + scale(cos(rad)) +
# scale(Day)+
scale(Temp) +
scale(StringencyIndex) +
(1|weekday) + (1 | genus) + (1 | Species) + (1 | sp_day_year) + (1 | IDLocality) + (1 | sp_loc),
# (1 | Year) + (1 | weekday) + (1|genus) + (1 | Species) + (1 | sp_day_year) + (1 | IDLocality),
data = s[Country == "Australia"], REML = FALSE
)
est_ass <- est_out(ass, "Australia: (1|weekday) + (1|genus)+(1|Species)+(1|sp_day_year)+(1|IDLocality)+(1|sp_loc)")
est_ass[, control_for_starting_distance := "yes"]
asx <- lmer(scale(log(FID)) ~
scale(Year) +
# scale(log(SD)) +
scale(log(FlockSize)) +
scale(log(BodyMass)) +
scale(sin(rad)) + scale(cos(rad)) +
# scale(Day)+
scale(Temp) +
scale(StringencyIndex) +
(1|weekday) + (1 | genus) + (1 | Species) + (1 | sp_day_year) + (1 | IDLocality) + (1 | sp_loc),
# (1 | Year) + (1 | weekday) + (1 | genus) + (1 | Species) + (1 | sp_day_year) + (1 | IDLocality),
data = s[Country == "Australia"], , REML = FALSE,
control = lmerControl(
optimizer ='optimx', optCtrl=list(method='nlminb'))
)
est_asx <- est_out(asx, "Australia: (1|weekday) + (1|genus)+(1|Species)+(1|sp_day_year)+(1|IDLocality)+(1|sp_loc)")
est_asx[, control_for_starting_distance := "no"]
# PL
pss <- lmer(scale(log(FID)) ~
scale(Year) +
scale(log(SD)) +
scale(log(FlockSize)) +
scale(log(BodyMass)) +
scale(sin(rad)) + scale(cos(rad)) +
# scale(Day)+
scale(Temp) +
scale(StringencyIndex) +
(1|weekday) + (1 | genus) + (1 | Species) + (1 | sp_day_year),
# (1 | Year) + (1 | weekday) + (1|genus) + (1 | Species) + (1 | sp_day_year) + (1 | IDLocality),
data = s[Country == "Poland"], REML = FALSE
)
est_pss <- est_out(pss, "Poland: (1|weekday) + (1|genus)+(1|Species)+(1|sp_day_year)")
est_pss[, control_for_starting_distance := "yes"]
psx <- lmer(scale(log(FID)) ~
scale(Year) +
# scale(log(SD)) +
scale(log(FlockSize)) +
scale(log(BodyMass)) +
scale(sin(rad)) + scale(cos(rad)) +
# scale(Day)+
scale(Temp) +
scale(StringencyIndex) +
(1|weekday) + (1 | genus) + (1 | Species) + (1 | sp_day_year),
# (1 | Year) + (1 | weekday) + (1 | genus) + (1 | Species) + (1 | sp_day_year) + (1 | IDLocality),
data = s[Country == "Poland"], REML = FALSE
)
est_psx <- est_out(psx, "Poland: (1|weekday) + (1|genus)+(1|Species)+(1|sp_day_year)")
est_psx[, control_for_starting_distance := "no"]
# combine
est_mss[, Country := 'All\n(mixed model)']
est_msx[, Country := "All\n(mixed model)"]
est_ass[, Country := "Australia"]
est_asx[, Country := "Australia"]
est_css[, Country := "Czechia"]
est_csx[, Country := "Czechia"]
est_hss[, Country := "Hungary"]
est_hsx[, Country := "Hungary"]
est_pss[, Country := "Poland"]
est_psx[, Country := "Poland"]
est_fss[, Country := "Finland"]
est_fsx[, Country := "Finland"]
os = rbind(est_mss, est_msx,
est_ass, est_asx,
est_css, est_csx,
est_hss, est_hsx,
est_pss, est_psx,
est_fss, est_fsx)
save(os, file = here::here('Data/dat_est_Stringency_rev.Rdata'))
# estimates for table
mss_out <- m_out(name = "Table S4 - full a", dep = "Escape distance", model = mss, nsim = 5000)
msx_out <- m_out(name = "Table S4 - full b", dep = "Escape distance", model = msx, nsim = 5000)
css_out <- m_out(name = "Table S4 - CZ a", dep = "Escape distance", model = css, nsim = 5000)
csx_out <- m_out(name = "Table S4 - CZ b", dep = "Escape distance", model = csx, nsim = 5000)
fss_out <- m_out(name = "Table S4 - FI a", dep = "Escape distance", model = fss, nsim = 5000)
fsx_out <- m_out(name = "Table S4 - FI b", dep = "Escape distance", model = fsx, nsim = 5000)
hss_out <- m_out(name = "Table S4 - HU a", dep = "Escape distance", model = hss, nsim = 5000)
hsx_out <- m_out(name = "Table S4 - HU b", dep = "Escape distance", model = hsx, nsim = 5000)
ass_out <- m_out(name = "Table S4 - AU a", dep = "Escape distance", model = ass, nsim = 5000)
asx_out <- m_out(name = "Table S4 - AU b", dep = "Escape distancey", model = asx, nsim = 5000)
pss_out <- m_out(name = "Table S4 - PL a", dep = "Escape distance", model = pss, nsim = 5000)
psx_out <- m_out(name = "Table S4 - PL b", dep = "Escape distancey", model = psx, nsim = 5000)
out_FID_s <- rbind(mss_out, msx_out, fss_out, fsx_out, pss_out, psx_out, css_out, csx_out, hss_out, hsx_out, ass_out, asx_out, fill = TRUE)
out_FID_s[is.na(out_FID_s)] <- ""
out_FID_s$R2_mar = out_FID_s$R2_con = NULL
out_FID_s[, effect := gsub("scale\\(StringencyIndex\\)", "stringency index", effect)]
out_FID_s[, effect := gsub("scale\\(Year\\)", "year", effect)]
out_FID_s[, effect := gsub("scale\\(log\\(SD\\)", "starting distance (ln)", effect)]
out_FID_s[, effect := gsub("scale\\(Temp\\)", "temperaturre", effect)]
out_FID_s[, effect := gsub("scale\\(log\\(FlockSize\\)\\)", "flock size (ln)", effect)]
out_FID_s[, effect := gsub("scale\\(log\\(BodyMass\\)\\)", "body mass (ln)", effect)]
out_FID_s[, effect := gsub("scale\\(sin\\(rad\\)\\)", "time (sine of radians)", effect)]
out_FID_s[, effect := gsub("scale\\(cos\\(rad\\)\\)", "time (cosine of radians)", effect)]
out_FID_s[, effect := gsub("Species", "species", effect)]
out_FID_s[, effect := gsub("Country", "country", effect)]
out_FID_s[, effect := gsub("sp_day_year", "species within day & year", effect)]
out_FID_s[, effect := gsub("IDLocality", "site", effect)]
out_FID_s[, effect := gsub("sp_loc", "species within site", effect)]
out_FID_s[type == "random" & grepl("stringency index", effect, fixed = TRUE), effect := paste("stringency index (slope) |", gsub(" stringency index", "", effect))]
fwrite(file = here::here("Outputs/Table_S4_rev.csv"), out_FID_s)
load(here::here("Data/dat_est_Stringency_rev.Rdata"))
os[predictor %in% c("scale(StringencyIndex)"), predictor := "Stringency Index"]
oso <- os[predictor %in% c("Stringency Index")]
oso[, N:=as.numeric(sub('.*N = ', '', model))]
# add meta-analytical mean
oso_s = oso[control_for_starting_distance == 'yes']
met = summary(meta.summaries(d = oso_s$estimate, se = oso_s$sd, method = "fixed", weights = oso_s$N))$summci
oso_met = data.table(predictor = "Period", estimate = met[2], lwr = met[1], upr = met[3], sd = NA, model = NA, control_for_starting_distance = "yes", Country = "Combined\n(metanalytical)", N = NA)
oso_sx = oso[control_for_starting_distance == "no"]
metx = summary(meta.summaries(d = oso_sx$estimate, se = oso_sx$sd, method = "fixed", weights = oso_sx$N))$summci
oso_metx = data.table(predictor = "Period", estimate = metx[2], lwr = metx[1], upr = metx[3], sd = NA, model = NA, control_for_starting_distance = "no", Country = "Combined\n(metanalytical)", N = NA)
oso = rbind(oso, oso_met, oso_metx)
oso[, Country := factor(Country, levels = rev(c("Finland", "Poland", "Czechia", "Hungary", "Australia", "Combined\n(metanalytical)", "All\n(mixed model)")))]
# prepare for adding N
oso[, N := as.character(N)]
oso[control_for_starting_distance == "no" | is.na(N), N := ""]
oso[, n_pos := 0.35]
width_ <- .5 # spacing between error bars
#col_ <- c(brewer.pal(n = 12, name = "Paired"), "grey30", "grey80")
#Tol_bright <- c("#EE6677", "#228833", "#4477AA", "#CCBB44", "#66CCEE", "#AA3377", "#BBBBBB")
#Tol_muted <- c("#88CCEE", "#44AA99", "#117733", "#332288", "#DDCC77", "#999933", "#CC6677", "#882255", "#AA4499", "#DDDDDD")
#Tol_light <- c("#BBCC33", "#AAAA00", "#77AADD", "#EE8866", "#EEDD88", "#FFAABB", "#99DDFF", "#44BB99", "#DDDDDD")
# From Color Universal Design (CUD): https://jfly.uni-koeln.de/color/
#Okabe_Ito <- c("#E69F00", "#56B4E9", "#009E73", "#F0E442", "#0072B2", "#D55E00", "#CC79A7", "#000000")
#col_ = Okabe_Ito[7:1]
# JAMA and LocusZoom modified order
#col_ = c("#374E55FF", "#374E55FF", "#DF8F44FF", "#79AF97FF", "#00A1D5FF", "#B24745FF", "#80796BFF") #"#6A6599FF",
#col_ <- c("#357EBDFF", "#9632B8FF", "#46B8DAFF", "#5CB85CFF", "#EEA236FF", "#D43F3AFF", "#D43F3AFF")[7:1] # "#D43F3AFF", "#B8B8B8FF"
col_ = c("#357EBDFF", "#D43F3AFF", "#46B8DAFF", "#5CB85CFF", "#EEA236FF", "#9632B8FF", "#9632B8FF")[7:1] # "#D43F3AFF", "#B8B8B8FF"
#show_col(col_)
gs6a =
ggplot(oso, aes(x = estimate, y = Country, col = Country, shape = control_for_starting_distance)) +
geom_vline(xintercept = 0, color = "grey", linetype = "dotted") +
geom_errorbarh(aes(xmin = lwr, xmax = upr), height = 0, position = ggstance::position_dodgev(width_)) +
# geom_point(position = ggstance::position_dodgev(.6)) +
geom_point(position = position_dodge(width = width_), bg = "white", size = 1.1) +
# scale_color_viridis(discrete=TRUE, begin=0, end = 0.5) +
# scale_fill_viridis(discrete=TRUE, begin=0, end = 0.5) +
# geom_text( aes(x = n_pos,label = N), vjust = 0, size = 1.75, position = ggstance::position_dodgev(width_))+ # 3 positions for 3 bars
# annotate("text", x=log10(3), y=85, label= "Used", col = "grey30", size = 2.5)+
geom_text( aes(x = n_pos,label = N), vjust = 1, size = 1.75, position = ggstance::position_dodgev(width_))+
scale_shape_manual(name = "Controlled for\nstarting distance", guide = guide_legend(reverse = TRUE), values = c(21, 19)) +
#scale_color_jama(guide = "none")+ #, palette = 'light'
scale_color_manual(guide = "none", values = col_) + #guide_legend(reverse = TRUE)
scale_x_continuous(breaks = round(seq(-0.3, 0.4, by = 0.1), 1)) +
ylab("") +
xlab("Standardized effect size of\nStringency Index\n[on flight initiation distance]") +
labs(tag = 'a)')+
# coord_cartesian(xlim = c(-.15, .15)) +
# scale_x_continuous(breaks = round(seq(-.15, .15, by = 0.05),2)) +
theme_bw() +
theme(
legend.position = "right",
plot.tag = element_text(size = 7),
legend.title = element_text(size = 7),
legend.text = element_text(size = 6),
# legend.spacing.y = unit(0.1, 'cm'),
legend.key.height = unit(0.5, "line"),
legend.margin = margin(0, 0, 0, 0),
# legend.position=c(0.5,1.6),
plot.title = element_text(color = "grey", size = 7),
plot.margin = margin(b = 0.5, l = 0.5, t = 0.5, r = 1, unit = "pt"),
panel.grid = element_blank(),
panel.border = element_blank(),
panel.background = element_blank(),
axis.line = element_line(colour = ax_lines, size = 0.25),
axis.line.y = element_blank(),
axis.ticks.y = element_blank(),
axis.ticks.x = element_line(colour = ax_lines, size = 0.25),
# axis.text.x = element_text()
axis.ticks.length = unit(1, "pt"),
axis.text.x = element_text(, size = 6),
axis.text.y = element_text(colour = "black", size = 7),
axis.title = element_text(size = 7)
)
if(save_plot==TRUE){
ggsave(here::here("Outputs/Fig_S6a_Stringency.png"), gs6a, width = 8, height = 6.5, unit = "cm", dpi = 600)
}
# predictions for Fig and Table - Google Mobility
# full
mgs <- lmer(scale(log(FID)) ~
scale(Year) +
scale(log(SD)) +
scale(log(FlockSize)) +
scale(log(BodyMass)) +
scale(sin(rad)) + scale(cos(rad)) +
# scale(Day)+
scale(Temp) +
scale(parks_percent_change_from_baseline) +
(1|weekday) + (1| genus) + (1 | Species) + (1 | sp_day_year) +
(scale(parks_percent_change_from_baseline)| Country) + (1 | IDLocality) + (1 | sp_loc),
data = ss, REML = FALSE,
control = lmerControl(
optimizer = "optimx", optCtrl = list(method = "nlminb")
)
)
est_mgs <- est_out(mgs, "ALL: (1|weekday) + (1|genus) + (1|Species) + (1|sp_day_year) + (scale(Google)|Country) + (1|IDLocality) +(1|sp_loc)")
est_mgs[, control_for_starting_distance := "yes"]
mgx <- lmer(scale(log(FID)) ~
scale(Year) +
#scale(log(SD)) +
scale(log(FlockSize)) +
scale(log(BodyMass)) +
scale(sin(rad)) + scale(cos(rad)) +
# scale(Day)+
scale(Temp) +
scale(parks_percent_change_from_baseline) +
(1|weekday) + (1| genus) + (1 | Species) + (1 | sp_day_year) +
(scale(parks_percent_change_from_baseline) | Country) + (1 | IDLocality) + (1 | sp_loc),
data = ss, REML = FALSE,
control = lmerControl(
optimizer = "optimx", optCtrl = list(method = "nlminb")
)
)
est_mgx <- est_out(mgx, "ALL: (1|weekday) + (1|genus) + (1|Species) + (1|sp_day_year) + (scale(Google)|Country) + (1|IDLocality) +(1|sp_loc)")
est_mgx[, control_for_starting_distance := "no"]
# CZ - singular fits only due to genera estimated as zero (removing it changes no results)
cgs <- lmer(scale(log(FID)) ~
#scale(Year) +
scale(log(SD)) +
scale(log(FlockSize)) +
scale(log(BodyMass)) +
scale(sin(rad)) + scale(cos(rad)) +
# scale(Day)+
scale(Temp) +
scale(parks_percent_change_from_baseline) +
(1|weekday) + (1 | genus) + (1| Species) + (1 | sp_day_year) + (1| IDLocality) + (1|sp_loc),
#(1 | Year) + (1 | weekday) + (1|genus) + (1 | Species) + (1 | sp_day_year) + (1 | IDLocality),
data = ss[Country == "Czechia"], REML = FALSE
)
est_cgs <- est_out(cgs, "Czechia:(1|weekday)+(1|genus)+(1|Species)+(1|sp_day_year)+(1|IDLocality)+(1|sp_loc)")
est_cgs[, control_for_starting_distance := "yes"]
cgx <- lmer(scale(log(FID)) ~
#scale(Year) +
scale(log(SD)) +
scale(log(FlockSize)) +
scale(log(BodyMass)) +
scale(sin(rad)) + scale(cos(rad)) +
# scale(Day)+
scale(Temp) +
scale(parks_percent_change_from_baseline) +
(1|weekday) + (1 | genus) + (1| Species) + (1 | sp_day_year) + (1| IDLocality) + (1|sp_loc),
#(1 | Year) + (1 | weekday) + (1|genus) + (1 | Species) + (1 | sp_day_year) + (1 | IDLocality),
data = ss[Country == "Czechia"], REML = FALSE
)
est_cgx <- est_out(cgx, "Czechia: (1|weekday) + (1|genus)+(1|Species)+(1|sp_day_year)+(1|IDLocality)+(1|sp_loc)")
est_cgx[, control_for_starting_distance := "no"]
# FI
fgs <- lmer(scale(log(FID)) ~
scale(Year) +
scale(log(SD)) +
scale(log(FlockSize)) +
scale(log(BodyMass)) +
scale(sin(rad)) + scale(cos(rad)) +
# scale(Day)+
scale(Temp) +
scale(parks_percent_change_from_baseline) +
(1|weekday) + (1 | genus) + (1| Species) + (1 | sp_day_year) + (1| IDLocality) + (1|sp_loc),
data = ss[Country == "Finland"], REML = FALSE
)
est_fgs <- est_out(fgs, "Finland: (1|weekday) + (1|genus)+(1|Species)+(1|sp_day_year)+(1|IDLocality)+(1|sp_loc)")
est_fgs[, control_for_starting_distance := "yes"]
fgx <- lmer(scale(log(FID)) ~
scale(Year) +
# scale(log(SD)) +
scale(log(FlockSize)) +
scale(log(BodyMass)) +
scale(sin(rad)) + scale(cos(rad)) +
# scale(Day)+
scale(Temp) +
scale(parks_percent_change_from_baseline) +
(1|weekday) + (1 | genus) + (1| Species) + (1 | sp_day_year) + (1| IDLocality) + (1|sp_loc),
data = ss[Country == "Finland"], REML = FALSE
)
est_fgx <- est_out(fgx, "Finland: (1|weekday) + (1|genus)+(1|Species)+(1|sp_day_year)+(1|IDLocality)+(1|sp_loc)")
est_fgx[, control_for_starting_distance := "no"]
# HU - singular fits only due to sp_loc estimated as zero (removing it changes no results)
hgs <- lmer(scale(log(FID)) ~
scale(Year) +
scale(log(SD)) +
scale(log(FlockSize)) +
scale(log(BodyMass)) +
scale(sin(rad)) + scale(cos(rad)) +
# scale(Day)+
scale(Temp) +
scale(parks_percent_change_from_baseline) +
(1|weekday) + (1 | genus) + (1 | Species) + (1 | sp_day_year) + (1 | IDLocality) + (1 | sp_loc),
# (1 | Year) + (1 | weekday) + (1|genus) + (1 | Species) + (1 | sp_day_year) + (1 | IDLocality),
data = ss[Country == "Hungary"], REML = FALSE
)
est_hgs <- est_out(hgs, "Hungary: (1|weekday) + (1|genus)+(1|Species)+(1|sp_day_year)+(1|IDLocality)+(1|sp_loc)")
est_hgs[, control_for_starting_distance := "yes"]
hgx <- lmer(scale(log(FID)) ~
scale(Year) +
# scale(log(SD)) +
scale(log(FlockSize)) +
scale(log(BodyMass)) +
scale(sin(rad)) + scale(cos(rad)) +
# scale(Day)+
scale(Temp) +
scale(parks_percent_change_from_baseline) +
(1|weekday) + (1 | genus) + (1 | Species) + (1 | sp_day_year) + (1 | IDLocality) + (1 | sp_loc),
# (1 | Year) + (1 | weekday) + (1 | genus) + (1 | Species) + (1 | sp_day_year) + (1 | IDLocality),
data = ss[Country == "Hungary"], REML = FALSE
)
est_hgx <- est_out(hgx, "Hungary: (1|weekday) + (1|genus)+(1|Species)+(1|sp_day_year)+(1|IDLocality)+(1|sp_loc)")
est_hgx[, control_for_starting_distance := "no"]
# AU - singular fits only due to Year and random slope estimated as zero (removing those changes no results)
ags <- lmer(scale(log(FID)) ~
scale(Year) +
scale(log(SD)) +
scale(log(FlockSize)) +
scale(log(BodyMass)) +
scale(sin(rad)) + scale(cos(rad)) +
# scale(Day)+
scale(Temp) +
scale(parks_percent_change_from_baseline) +
(1|weekday) + (1 | genus) + (1 | Species) + (1 | sp_day_year) + (1 | IDLocality) + (1 | sp_loc),
# (1 | Year) + (1 | weekday) + (1|genus) + (1 | Species) + (1 | sp_day_year) + (1 | IDLocality),
data = ss[Country == "Australia"], REML = FALSE
)
est_ags <- est_out(ags, "Australia: (1|weekday) + (1|genus)+(1|Species)+(1|sp_day_year)+(1|IDLocality)+(1|sp_loc)")
est_ags[, control_for_starting_distance := "yes"]
agx <- lmer(scale(log(FID)) ~
scale(Year) +
# scale(log(SD)) +
scale(log(FlockSize)) +
scale(log(BodyMass)) +
scale(sin(rad)) + scale(cos(rad)) +
# scale(Day)+
scale(Temp) +
scale(parks_percent_change_from_baseline) +
(1|weekday) + (1 | genus) + (1 | Species) + (1 | sp_day_year) + (1 | IDLocality) + (1 | sp_loc),
# (1 | Year) + (1 | weekday) + (1 | genus) + (1 | Species) + (1 | sp_day_year) + (1 | IDLocality),
data = ss[Country == "Australia"], , REML = FALSE,
control = lmerControl(
optimizer ='optimx', optCtrl=list(method='nlminb'))
)
est_agx <- est_out(agx, "Australia: (1|weekday) + (1|genus)+(1|Species)+(1|sp_day_year)+(1|IDLocality)+(1|sp_loc)")
est_agx[, control_for_starting_distance := "no"]
# PL
pgs <- lmer(scale(log(FID)) ~
scale(Year) +
scale(log(SD)) +
scale(log(FlockSize)) +
scale(log(BodyMass)) +
scale(sin(rad)) + scale(cos(rad)) +
# scale(Day)+
scale(Temp) +
scale(parks_percent_change_from_baseline) +
(1|weekday) + (1 | genus) + (1 | Species) + (1 | sp_day_year),
# (1 | Year) + (1 | weekday) + (1|genus) + (1 | Species) + (1 | sp_day_year) + (1 | IDLocality),
data = ss[Country == "Poland"], REML = FALSE
)
est_pgs <- est_out(pgs, "Poland: (1|weekday) + (1|genus)+(1|Species)+(1|sp_day_year)")
est_pgs[, control_for_starting_distance := "yes"]
pgx <- lmer(scale(log(FID)) ~
scale(Year) +
# scale(log(SD)) +
scale(log(FlockSize)) +
scale(log(BodyMass)) +
scale(sin(rad)) + scale(cos(rad)) +
# scale(Day)+
scale(Temp) +
scale(parks_percent_change_from_baseline) +
(1|weekday) + (1 | genus) + (1 | Species)+ (1 | sp_day_year),
# (1 | Year) + (1 | weekday) + (1 | genus) + (1 | Species) + (1 | sp_day_year) + (1 | IDLocality),
data = ss[Country == "Poland"], REML = FALSE
)
est_pgx <- est_out(pgx, "Poland: (1|weekday) + (1|genus)+(1|Species)+(1|sp_day_year)")
est_pgx[, control_for_starting_distance := "no"]
# combine
est_mgs[, Country := 'All\n(mixed model)']
est_mgx[, Country := "All\n(mixed model)"]
est_ags[, Country := "Australia"]
est_agx[, Country := "Australia"]
est_cgs[, Country := "Czechia"]
est_cgx[, Country := "Czechia"]
est_hgs[, Country := "Hungary"]
est_hgx[, Country := "Hungary"]
est_pgs[, Country := "Poland"]
est_pgx[, Country := "Poland"]
est_fgs[, Country := "Finland"]
est_fgx[, Country := "Finland"]
og = rbind(est_mgs, est_mgx,
est_ags, est_agx,
est_cgs, est_cgx,
est_hgs, est_hgx,
est_pgs, est_pgx,
est_fgs, est_fgx)
save(og, file = here::here('Data/dat_est_Google_rev.Rdata'))
# estimatees for table
mgs_out <- m_out(name = "Table S5 - full a", dep = "Escape distance", model = mgs, nsim = 5000)
mgx_out <- m_out(name = "Table S5 - full b", dep = "Escape distance", model = mgx, nsim = 5000)
cgs_out <- m_out(name = "Table S5 - CZ a", dep = "Escape distance", model = cgs, nsim = 5000)
cgx_out <- m_out(name = "Table S5 - CZ b", dep = "Escape distance", model = cgx, nsim = 5000)
fgs_out <- m_out(name = "Table S5 - FI a", dep = "Escape distance", model = fgs, nsim = 5000)
fgx_out <- m_out(name = "Table S5 - FI b", dep = "Escape distance", model = fgx, nsim = 5000)
hgs_out <- m_out(name = "Table S5 - HU a", dep = "Escape distance", model = hgs, nsim = 5000)
hgx_out <- m_out(name = "Table S5 - HU b", dep = "Escape distance", model = hgx, nsim = 5000)
ags_out <- m_out(name = "Table S5 - AU a", dep = "Escape distance", model = ags, nsim = 5000)
agx_out <- m_out(name = "Table S5 - AU b", dep = "Escape distancey", model = agx, nsim = 5000)
pgs_out <- m_out(name = "Table S5 - PL a", dep = "Escape distance", model = pgs, nsim = 5000)
pgx_out <- m_out(name = "Table S5 - PL b", dep = "Escape distancey", model = pgx, nsim = 5000)
out_FID_g <- rbind(mgs_out, mgx_out, fgs_out, fgx_out, pgs_out, pgx_out, cgs_out, cgx_out, hgs_out, hgx_out, ags_out, agx_out, fill = TRUE)
out_FID_g[is.na(out_FID_g)] <- ""
out_FID_g[, effect := gsub("scale\\(parks_percent_change_from_baseline\\)", "Google Mobility", effect)]
out_FID_g[, effect := gsub("scale\\(Year\\)", "year", effect)]
out_FID_g[, effect := gsub("scale\\(log\\(SD\\)", "starting distance (ln)", effect)]
out_FID_g[, effect := gsub("scale\\(Temp\\)", "temperaturre", effect)]
out_FID_g[, effect := gsub("scale\\(log\\(FlockSize\\)\\)", "flock size (ln)", effect)]
out_FID_g[, effect := gsub("scale\\(log\\(BodyMass\\)\\)", "body mass (ln)", effect)]
out_FID_g[, effect := gsub("scale\\(sin\\(rad\\)\\)", "time (sine of radians)", effect)]
out_FID_g[, effect := gsub("scale\\(cos\\(rad\\)\\)", "time (cosine of radians)", effect)]
out_FID_g[, effect := gsub("Species", "species", effect)]
out_FID_g[, effect := gsub("Country", "country", effect)]
out_FID_g[, effect := gsub("sp_day_year", "species within day & year", effect)]
out_FID_g[, effect := gsub("IDLocality", "site", effect)]
out_FID_g[, effect := gsub("sp_loc", "species within site", effect)]
out_FID_g[, effect := paste("Google Mobility (slope) |", gsub(" Google Mobility", "", effect))]
fwrite(file = here::here("Outputs/Table_S5_rev.csv"), out_FID_g)
load(here::here("Data/dat_est_Google_rev.Rdata"))
og[predictor %in% c("scale(parks_percent_change_from_baseline)"), predictor := "Google Mobility"]
ogo <- og[predictor %in% c("Google Mobility")]
ogo[, N:=as.numeric(sub('.*N = ', '', model))]
# add meta-analytical mean
ogo_s = ogo[control_for_starting_distance == 'yes']
met = summary(meta.summaries(d = ogo_s$estimate, se = ogo_s$sd, method = "fixed", weights = ogo_s$N))$summci
ogo_met = data.table(predictor = "Period", estimate = met[2], lwr = met[1], upr = met[3], sd = NA, model = NA, control_for_starting_distance = "yes", Country = "Combined\n(metanalytical)", N = NA)
ogo_sx = ogo[control_for_starting_distance == "no"]
metx = summary(meta.summaries(d = ogo_sx$estimate, se = ogo_sx$sd, method = "fixed", weights = ogo_sx$N))$summci
ogo_metx = data.table(predictor = "Period", estimate = metx[2], lwr = metx[1], upr = metx[3], sd = NA, model = NA, control_for_starting_distance = "no", Country = "Combined\n(metanalytical)", N = NA)
ogo = rbind(ogo, ogo_met, ogo_metx)
ogo[, Country := factor(Country, levels = rev(c("Finland", "Poland", "Czechia", "Hungary", "Australia", "Combined\n(metanalytical)", "All\n(mixed model)")))]
# prepare for adding N
ogo[, N := as.character(N)]
ogo[control_for_starting_distance == "no" | is.na(N), N := ""]
ogo[, n_pos := .15]
width_ <- .5 # spacing between error bars
#col_ <- c(brewer.pal(n = 12, name = "Paired"), "grey30", "grey80")
#Tol_bright <- c("#EE6677", "#228833", "#4477AA", "#CCBB44", "#66CCEE", "#AA3377", "#BBBBBB")
#Tol_muted <- c("#88CCEE", "#44AA99", "#117733", "#332288", "#DDCC77", "#999933", "#CC6677", "#882255", "#AA4499", "#DDDDDD")
#Tol_light <- c("#BBCC33", "#AAAA00", "#77AADD", "#EE8866", "#EEDD88", "#FFAABB", "#99DDFF", "#44BB99", "#DDDDDD")
# From Color Universal Design (CUD): https://jfly.uni-koeln.de/color/
#Okabe_Ito <- c("#E69F00", "#56B4E9", "#009E73", "#F0E442", "#0072B2", "#D55E00", "#CC79A7", "#000000")
#col_ = Okabe_Ito[7:1]
# JAMA and LocusZoom modified order
#col_ = c("#374E55FF", "#374E55FF", "#DF8F44FF", "#79AF97FF", "#00A1D5FF", "#B24745FF", "#80796BFF") #"#6A6599FF",
#col_ <- c("#357EBDFF", "#9632B8FF", "#46B8DAFF", "#5CB85CFF", "#EEA236FF", "#D43F3AFF", "#D43F3AFF")[7:1] # "#D43F3AFF", "#B8B8B8FF"
col_ = c("#357EBDFF", "#D43F3AFF", "#46B8DAFF", "#5CB85CFF", "#EEA236FF", "#9632B8FF", "#9632B8FF")[7:1] # "#D43F3AFF", "#B8B8B8FF"
#show_col(col_)
gs6b =
ggplot(ogo, aes(x = estimate, y = Country, col = Country, shape = control_for_starting_distance)) +
geom_vline(xintercept = 0, color = "grey", linetype = "dotted") +
geom_errorbarh(aes(xmin = lwr, xmax = upr), height = 0, position = ggstance::position_dodgev(width_)) +
# geom_point(position = ggstance::position_dodgev(.6)) +
geom_point(position = position_dodge(width = width_), bg = "white", size = 1.1) +
# scale_color_viridis(discrete=TRUE, begin=0, end = 0.5) +
# scale_fill_viridis(discrete=TRUE, begin=0, end = 0.5) +
# geom_text( aes(x = n_pos,label = N), vjust = 0, size = 1.75, position = ggstance::position_dodgev(width_))+ # 3 positions for 3 bars
# annotate("text", x=log10(3), y=85, label= "Used", col = "grey30", size = 2.5)+
geom_text( aes(x = n_pos,label = N), vjust = 1, size = 1.75, position = ggstance::position_dodgev(width_))+
scale_shape_manual(name = "Controlled for\nstarting distance", guide = guide_legend(reverse = TRUE), values = c(21, 19)) +
#scale_color_jama(guide = "none")+ #, palette = 'light'
scale_color_manual(guide = "none", values = col_) + #guide_legend(reverse = TRUE)
scale_x_continuous(breaks = round(seq(-0.3, 0.2, by = 0.1), 1)) +
ylab("") +
xlab("Standardized effect size of\nGoogle Mobility (human presence)\n[on flight initiation distance]") +
labs(tag = 'b)')+
# coord_cartesian(xlim = c(-.15, .15)) +
# scale_x_continuous(breaks = round(seq(-.15, .15, by = 0.05),2)) +
theme_bw() +
theme(
plot.tag = element_text(size = 7),
legend.position = "right",
legend.title = element_text(size = 7),
legend.text = element_text(size = 6),
# legend.spacing.y = unit(0.1, 'cm'),
legend.key.height = unit(0.5, "line"),
legend.margin = margin(0, 0, 0, 0),
# legend.position=c(0.5,1.6),
plot.title = element_text(color = "grey", size = 7),
plot.margin = margin(b = 0.5, l = 0.5, t = 0.5, r = 1, unit = "pt"),
panel.grid = element_blank(),
panel.border = element_blank(),
panel.background = element_blank(),
axis.line = element_line(colour = ax_lines, size = 0.25),
axis.line.y = element_blank(),
axis.ticks.y = element_blank(),
axis.ticks.x = element_line(colour = ax_lines, size = 0.25),
# axis.text.x = element_text()
axis.ticks.length = unit(1, "pt"),
axis.text.x = element_text(, size = 6),
axis.text.y = element_text(colour = "black", size = 7),
axis.title = element_text(size = 7)
)
if(save_plot==TRUE){
ggsave(here::here("Outputs/Fig_S6b_Google_rev_width_CustomLocusZoom_v2.png"), gs6b, width = 8, height = 6.5, unit = "cm", dpi = 600)
}#gg_S6 <- ggarrange(
# gs6a+theme(legend.position = "none"), gs6b,
#nrow = 1, widths=c(1, 1.3)
#)
gg_S6 <- ggarrange(
gs6a, NULL, gs6b, widths = c(1,0.05,1),
nrow = 1, common.legend = TRUE, legend = 'right'
)
if (save_plot == TRUE) {
ggsave(here::here("Outputs/Fig_S6_rev.png"), gg_S6, width = 8*1.5, height = 6.5, unit = "cm", dpi = 600)
}
gg_S6
Figure S6 | Changes in avian tolerance towards humans in response to (a) stringency of governmental measures, and (b) Google Mobility. The dots with horizontal lines represent estimated standardised effect size and their 95% confidence intervals, the numbers sample sizes. For the country-specific and “All”, the effect sizes and 95% confidence intervals come from the joint posterior distribution of 5000 simulated values generated by the sim function from the arm package (Gelman et al. 2022) using the mixed model outputs controlled for starting distance of the observer (filled circles) or not (empty circles; Table S4 and S5). The models were further controlled for year, flock size (ln-transfomred), body size (ln-transformed), temperature (also a proxy for a day within the breeding season: rPearson = 0.48; Fig. S2), and time of a day, as well as for the non-independence of data points by fitting random intercepts of weekday, genus, species, species at a given day and year, country (in All - a global mixed model), site, and species within a site, while fitting, in case of “All” (global model an all countries) stringency index or Google Mobility as a random slope within country (i.e. allowing for different effect in each country). Fitting stringency index or Google Mobility as random slope at other random intercepts produces similar results (Fig. S1b and S1c, Table S2b and S2c). The multicollinearity was small as correlations between predictors were weak (Fig. S2). For the “Combined”, the estimate and 95% confidence interval represent the meta-analytical mean based on the country estimates and their standard deviation (from the country-specific models), and sample size per country. These analyses were restricted to data collected in the Period during the COVID-19 shutdowns. Stringency index data were sourced from (Hale et al. 2021), Google Mobility from Google Mobility Reports.
Table S4 | Escape distance in relations to stringency index
out_FID_s$error_structure = out_FID_s$response = NULL
out_FID_s[model != "", model := c(
"All countries", "All countries, without starting distance",
"Finland", "Finland, without starting distance",
"Poland", "Poland, without starting distance",
"Czechia", "Czechia, without starting distance",
"Hungary", "Hungary, without starting distance",
"Australia", "Australia, without starting distance"
)]
setnames(out_FID_s, old = c("estimate_r", "lwr_r", "upr_r"), new = c("estimate", "lower", "upper"))
out_FID_s %>%
kbl() %>%
kable_paper("hover", full_width = F)| model | N | type | effect | estimate | lower | upper |
|---|---|---|---|---|---|---|
| All countries | 3676 | fixed | (Intercept) | -0.079 | -0.338 | 0.18 |
| fixed | year | 0.019 | -0.029 | 0.067 | ||
| fixed | starting distance (ln)) | 0.491 | 0.458 | 0.523 | ||
| fixed | flock size (ln) | -0.003 | -0.029 | 0.023 | ||
| fixed | body mass (ln) | 0.023 | -0.044 | 0.092 | ||
| fixed | time (sine of radians) | -0.02 | -0.061 | 0.02 | ||
| fixed | time (cosine of radians) | -0.005 | -0.041 | 0.03 | ||
| fixed | temperaturre | -0.007 | -0.048 | 0.035 | ||
| fixed | stringency index | 0.018 | -0.152 | 0.181 | ||
| random | species within day & year (Intercept) | 5% | 3% | 7% | ||
| random | species within site (Intercept) | 5% | 3% | 6% | ||
| random | species (Intercept) | 7% | 6% | 9% | ||
| random | site (Intercept) | 8% | 7% | 9% | ||
| random | genus (Intercept) | 4% | 4% | 4% | ||
| random | weekday (Intercept) | 0% | 0% | 0% | ||
| random | country (Intercept) | 2% | -16% | 16% | ||
| random | stringency index (slope) | country | 2% | -16% | 16% | ||
| random | Residual | 67% | 45% | 96% | ||
| All countries, without starting distance | 3676 | fixed | (Intercept) | 0.052 | -0.096 | 0.202 |
| fixed | year | 0.084 | 0.03 | 0.139 | ||
| fixed | flock size (ln) | 0.011 | -0.019 | 0.04 | ||
| fixed | body mass (ln) | 0.14 | 0.055 | 0.226 | ||
| fixed | time (sine of radians) | 0.046 | 0 | 0.093 | ||
| fixed | time (cosine of radians) | 0.049 | 0.009 | 0.091 | ||
| fixed | temperaturre | 0.005 | -0.042 | 0.051 | ||
| fixed | stringency index | -0.002 | -0.224 | 0.223 | ||
| random | species within day & year (Intercept) | 6% | 5% | 8% | ||
| random | species within site (Intercept) | 5% | 4% | 6% | ||
| random | species (Intercept) | 9% | 9% | 10% | ||
| random | site (Intercept) | 13% | 13% | 14% | ||
| random | genus (Intercept) | 5% | 5% | 6% | ||
| random | weekday (Intercept) | 0% | 0% | 0% | ||
| random | country (Intercept) | 3% | -9% | 9% | ||
| random | stringency index (slope) | country | 3% | -9% | 9% | ||
| random | Residual | 55% | 45% | 73% | ||
| Finland | 354 | fixed | (Intercept) | 0.032 | -0.186 | 0.237 |
| fixed | year | -0.108 | -0.342 | 0.119 | ||
| fixed | starting distance (ln)) | 0.189 | 0.088 | 0.288 | ||
| fixed | flock size (ln) | -0.111 | -0.211 | -0.011 | ||
| fixed | body mass (ln) | 0.175 | 0.017 | 0.333 | ||
| fixed | time (sine of radians) | 0.059 | -0.061 | 0.187 | ||
| fixed | time (cosine of radians) | 0.013 | -0.11 | 0.134 | ||
| fixed | temperaturre | -0.085 | -0.2 | 0.027 | ||
| fixed | stringency index | 0.154 | -0.07 | 0.375 | ||
| random | species within site (Intercept) | 12% | 12% | 13% | ||
| random | species within day & year (Intercept) | 3% | 3% | 3% | ||
| random | site (Intercept) | 12% | 10% | 15% | ||
| random | species (Intercept) | 0% | 0% | 1% | ||
| random | genus (Intercept) | 3% | 2% | 5% | ||
| random | weekday (Intercept) | 0% | 0% | 0% | ||
| random | Residual | 68% | 64% | 73% | ||
| Finland, without starting distance | 354 | fixed | (Intercept) | -0.008 | -0.237 | 0.217 |
| fixed | year | -0.066 | -0.297 | 0.166 | ||
| fixed | flock size (ln) | -0.103 | -0.205 | -0.002 | ||
| fixed | body mass (ln) | 0.245 | 0.085 | 0.409 | ||
| fixed | time (sine of radians) | 0.054 | -0.073 | 0.183 | ||
| fixed | time (cosine of radians) | 0.006 | -0.123 | 0.129 | ||
| fixed | temperaturre | -0.088 | -0.2 | 0.023 | ||
| fixed | stringency index | 0.135 | -0.088 | 0.36 | ||
| random | species within site (Intercept) | 14% | 14% | 15% | ||
| random | species within day & year (Intercept) | 4% | 3% | 4% | ||
| random | site (Intercept) | 14% | 11% | 16% | ||
| random | species (Intercept) | 3% | 2% | 5% | ||
| random | genus (Intercept) | 1% | 1% | 2% | ||
| random | weekday (Intercept) | 0% | 0% | 0% | ||
| random | Residual | 63% | 59% | 69% | ||
| Poland | 329 | fixed | (Intercept) | 0.026 | -0.197 | 0.254 |
| fixed | year | -0.206 | -0.476 | 0.068 | ||
| fixed | starting distance (ln)) | 0.533 | 0.447 | 0.619 | ||
| fixed | flock size (ln) | -0.032 | -0.113 | 0.049 | ||
| fixed | body mass (ln) | -0.083 | -0.186 | 0.027 | ||
| fixed | time (sine of radians) | 0.123 | -0.137 | 0.389 | ||
| fixed | time (cosine of radians) | -0.238 | -0.496 | 0.02 | ||
| fixed | temperaturre | -0.192 | -0.312 | -0.074 | ||
| fixed | stringency index | 0.041 | -0.228 | 0.31 | ||
| random | species within day & year (Intercept) | 7% | 6% | 7% | ||
| random | species (Intercept) | 9% | 7% | 10% | ||
| random | genus (Intercept) | 0% | 0% | 0% | ||
| random | weekday (Intercept) | 10% | 5% | 16% | ||
| random | Residual | 75% | 68% | 82% | ||
| Poland, without starting distance | 329 | fixed | (Intercept) | 0.023 | -0.242 | 0.285 |
| fixed | year | -0.341 | -0.65 | -0.023 | ||
| fixed | flock size (ln) | 0.053 | -0.045 | 0.148 | ||
| fixed | body mass (ln) | -0.011 | -0.163 | 0.145 | ||
| fixed | time (sine of radians) | -0.033 | -0.356 | 0.277 | ||
| fixed | time (cosine of radians) | -0.078 | -0.389 | 0.245 | ||
| fixed | temperaturre | -0.185 | -0.321 | -0.043 | ||
| fixed | stringency index | 0.033 | -0.28 | 0.346 | ||
| random | species within day & year (Intercept) | 4% | 4% | 4% | ||
| random | species (Intercept) | 21% | 17% | 25% | ||
| random | genus (Intercept) | 0% | 0% | 0% | ||
| random | weekday (Intercept) | 5% | 2% | 9% | ||
| random | Residual | 69% | 62% | 77% | ||
| Czechia | 1054 | fixed | (Intercept) | 0.213 | -0.008 | 0.437 |
| fixed | starting distance (ln)) | 0.344 | 0.283 | 0.405 | ||
| fixed | flock size (ln) | 0.005 | -0.045 | 0.055 | ||
| fixed | body mass (ln) | 0.183 | 0.025 | 0.34 | ||
| fixed | time (sine of radians) | 0.066 | -0.009 | 0.14 | ||
| fixed | time (cosine of radians) | 0.003 | -0.069 | 0.076 | ||
| fixed | temperaturre | 0.022 | -0.073 | 0.116 | ||
| fixed | stringency index | -0.052 | -0.173 | 0.069 | ||
| random | species within site (Intercept) | 4% | 4% | 4% | ||
| random | species within day & year (Intercept) | 2% | 2% | 2% | ||
| random | species (Intercept) | 14% | 12% | 17% | ||
| random | genus (Intercept) | 14% | 10% | 17% | ||
| random | site (Intercept) | 5% | 4% | 6% | ||
| random | weekday (Intercept) | 1% | 0% | 2% | ||
| random | Residual | 61% | 53% | 69% | ||
| Czechia, without starting distance | 1054 | fixed | (Intercept) | 0.216 | -0.006 | 0.435 |
| fixed | starting distance (ln)) | 0.345 | 0.286 | 0.405 | ||
| fixed | flock size (ln) | 0.006 | -0.044 | 0.055 | ||
| fixed | body mass (ln) | 0.185 | 0.028 | 0.335 | ||
| fixed | time (sine of radians) | 0.064 | -0.01 | 0.14 | ||
| fixed | time (cosine of radians) | 0.001 | -0.07 | 0.076 | ||
| fixed | temperaturre | 0.024 | -0.071 | 0.119 | ||
| fixed | stringency index | -0.052 | -0.172 | 0.071 | ||
| random | species within site (Intercept) | 4% | 4% | 4% | ||
| random | species within day & year (Intercept) | 2% | 2% | 2% | ||
| random | species (Intercept) | 14% | 12% | 17% | ||
| random | genus (Intercept) | 14% | 10% | 17% | ||
| random | site (Intercept) | 5% | 4% | 6% | ||
| random | weekday (Intercept) | 1% | 0% | 2% | ||
| random | Residual | 61% | 53% | 69% | ||
| Hungary | 874 | fixed | (Intercept) | 0.065 | -0.222 | 0.35 |
| fixed | year | 0.165 | 0.049 | 0.275 | ||
| fixed | starting distance (ln)) | 0.505 | 0.442 | 0.565 | ||
| fixed | flock size (ln) | -0.001 | -0.054 | 0.052 | ||
| fixed | body mass (ln) | 0.003 | -0.142 | 0.149 | ||
| fixed | time (sine of radians) | -0.126 | -0.195 | -0.057 | ||
| fixed | time (cosine of radians) | 0.024 | -0.051 | 0.101 | ||
| fixed | temperaturre | -0.018 | -0.124 | 0.082 | ||
| fixed | stringency index | 0.031 | -0.121 | 0.179 | ||
| random | species within day & year (Intercept) | 2% | 2% | 2% | ||
| random | species within site (Intercept) | 0% | 0% | 0% | ||
| random | species (Intercept) | 1% | 1% | 1% | ||
| random | genus (Intercept) | 12% | 8% | 16% | ||
| random | weekday (Intercept) | 1% | 0% | 2% | ||
| random | site (Intercept) | 10% | 6% | 13% | ||
| random | Residual | 74% | 65% | 83% | ||
| Hungary, without starting distance | 874 | fixed | (Intercept) | 0.06 | -0.431 | 0.552 |
| fixed | year | 0.205 | 0.073 | 0.334 | ||
| fixed | flock size (ln) | 0.002 | -0.058 | 0.063 | ||
| fixed | body mass (ln) | 0.299 | 0.067 | 0.521 | ||
| fixed | time (sine of radians) | -0.148 | -0.226 | -0.069 | ||
| fixed | time (cosine of radians) | 0.028 | -0.058 | 0.114 | ||
| fixed | temperaturre | -0.051 | -0.173 | 0.071 | ||
| fixed | stringency index | -0.005 | -0.175 | 0.168 | ||
| random | species within day & year (Intercept) | 3% | 3% | 4% | ||
| random | species within site (Intercept) | 0% | 0% | 0% | ||
| random | species (Intercept) | 1% | 1% | 1% | ||
| random | genus (Intercept) | 27% | 20% | 32% | ||
| random | weekday (Intercept) | 1% | 0% | 2% | ||
| random | site (Intercept) | 17% | 14% | 19% | ||
| random | Residual | 51% | 42% | 61% | ||
| Australia | 1065 | fixed | (Intercept) | -0.05 | -0.187 | 0.092 |
| fixed | year | 0.037 | -0.027 | 0.101 | ||
| fixed | starting distance (ln)) | 0.496 | 0.442 | 0.549 | ||
| fixed | flock size (ln) | 0.024 | -0.024 | 0.074 | ||
| fixed | body mass (ln) | -0.05 | -0.148 | 0.042 | ||
| fixed | time (sine of radians) | -0.037 | -0.11 | 0.036 | ||
| fixed | time (cosine of radians) | -0.051 | -0.118 | 0.015 | ||
| fixed | temperaturre | 0.033 | -0.027 | 0.097 | ||
| fixed | stringency index | 0.058 | -0.009 | 0.12 | ||
| random | species within day & year (Intercept) | 6% | 6% | 6% | ||
| random | species within site (Intercept) | 12% | 12% | 12% | ||
| random | species (Intercept) | 11% | 9% | 12% | ||
| random | genus (Intercept) | 2% | 2% | 3% | ||
| random | site (Intercept) | 7% | 5% | 9% | ||
| random | weekday (Intercept) | 0% | 0% | 0% | ||
| random | Residual | 62% | 58% | 66% | ||
| Australia, without starting distance | 1065 | fixed | (Intercept) | -0.102 | -0.285 | 0.083 |
| fixed | year | 0.101 | 0.028 | 0.176 | ||
| fixed | flock size (ln) | 0.055 | -0.001 | 0.112 | ||
| fixed | body mass (ln) | 0.063 | -0.065 | 0.185 | ||
| fixed | time (sine of radians) | -0.01 | -0.096 | 0.072 | ||
| fixed | time (cosine of radians) | -0.021 | -0.099 | 0.055 | ||
| fixed | temperaturre | 0.04 | -0.032 | 0.112 | ||
| fixed | stringency index | 0.09 | 0.014 | 0.166 | ||
| random | species within day & year (Intercept) | 7% | 7% | 8% | ||
| random | species within site (Intercept) | 11% | 11% | 11% | ||
| random | species (Intercept) | 13% | 11% | 14% | ||
| random | genus (Intercept) | 7% | 6% | 8% | ||
| random | site (Intercept) | 8% | 6% | 10% | ||
| random | weekday (Intercept) | 0% | 0% | 0% | ||
| random | Residual | 54% | 49% | 59% |
Table S5 | Escape distance in relations to Google Mobilitty
out_FID_g$R2_mar = out_FID_g$R2_con = out_FID_g$error_structure = out_FID_g$response = NULL
out_FID_g[model != "", model := c(
"All countries", "All countries, without starting distance",
"Finland", "Finland, without starting distance",
"Poland", "Poland, without starting distance",
"Czechia", "Czechia, without starting distance",
"Hungary", "Hungary, without starting distance",
"Australia", "Australia, without starting distance"
)]
setnames(out_FID_g, old = c('estimate_r','lwr_r','upr_r'), new =c('estimate','lower','upper' ))
out_FID_g %>%
kbl() %>%
kable_paper("hover", full_width = F)| model | N | type | effect | estimate | lower | upper |
|---|---|---|---|---|---|---|
| All countries | 3644 | fixed | (Intercept) | 0.023 | -0.098 | 0.141 |
| fixed | year | 0.018 | -0.023 | 0.06 | ||
| fixed | starting distance (ln)) | 0.498 | 0.467 | 0.529 | ||
| fixed | flock size (ln) | -0.001 | -0.028 | 0.025 | ||
| fixed | body mass (ln) | 0.025 | -0.044 | 0.091 | ||
| fixed | time (sine of radians) | -0.009 | -0.049 | 0.033 | ||
| fixed | time (cosine of radians) | 0.011 | -0.023 | 0.045 | ||
| fixed | temperaturre | 0.021 | -0.016 | 0.058 | ||
| fixed | Google Mobility | -0.101 | -0.226 | 0.025 | ||
| random | species within day & year (Intercept) | 6% | 5% | 6% | ||
| random | species within site (Intercept) | 5% | 4% | 5% | ||
| random | species (Intercept) | 7% | 7% | 8% | ||
| random | site (Intercept) | 9% | 8% | 10% | ||
| random | genus (Intercept) | 3% | 3% | 4% | ||
| random | weekday (Intercept) | 0% | 0% | 0% | ||
| random | country (Intercept) | 1% | -3% | 4% | ||
| random | country Google Mobility | 1% | -3% | 4% | ||
| random | Residual | 68% | 61% | 78% | ||
| All countries, without starting distance | 3644 | fixed | (Intercept) | 0.1 | -0.114 | 0.317 |
| fixed | year | 0.071 | 0.02 | 0.119 | ||
| fixed | flock size (ln) | 0.012 | -0.018 | 0.042 | ||
| fixed | body mass (ln) | 0.131 | 0.046 | 0.221 | ||
| fixed | time (sine of radians) | 0.062 | 0.015 | 0.109 | ||
| fixed | time (cosine of radians) | 0.074 | 0.037 | 0.112 | ||
| fixed | temperaturre | 0.038 | -0.005 | 0.083 | ||
| fixed | Google Mobility | -0.121 | -0.304 | 0.062 | ||
| random | species within day & year (Intercept) | 7% | 5% | 9% | ||
| random | species within site (Intercept) | 6% | 5% | 7% | ||
| random | species (Intercept) | 11% | 11% | 12% | ||
| random | site (Intercept) | 12% | 11% | 13% | ||
| random | genus (Intercept) | 5% | 4% | 5% | ||
| random | weekday (Intercept) | 0% | 0% | 0% | ||
| random | country (Intercept) | 1% | -13% | 9% | ||
| random | country Google Mobility | 1% | -13% | 9% | ||
| random | Residual | 59% | 46% | 79% | ||
| Finland | 322 | fixed | (Intercept) | 0.022 | -0.177 | 0.222 |
| fixed | year | 0.13 | -0.009 | 0.263 | ||
| fixed | starting distance (ln)) | 0.194 | 0.09 | 0.299 | ||
| fixed | flock size (ln) | -0.119 | -0.222 | -0.017 | ||
| fixed | body mass (ln) | 0.241 | 0.08 | 0.398 | ||
| fixed | time (sine of radians) | 0.132 | -0.013 | 0.276 | ||
| fixed | time (cosine of radians) | 0.044 | -0.081 | 0.166 | ||
| fixed | temperaturre | -0.005 | -0.136 | 0.127 | ||
| fixed | Google Mobility | -0.099 | -0.242 | 0.042 | ||
| random | species within site (Intercept) | 13% | 12% | 13% | ||
| random | species within day & year (Intercept) | 3% | 2% | 3% | ||
| random | site (Intercept) | 16% | 13% | 19% | ||
| random | species (Intercept) | 0% | 0% | 0% | ||
| random | genus (Intercept) | 2% | 1% | 4% | ||
| random | weekday (Intercept) | 0% | 0% | 0% | ||
| random | Residual | 66% | 62% | 71% | ||
| Finland, without starting distance | 322 | fixed | (Intercept) | -0.005 | -0.232 | 0.214 |
| fixed | year | 0.149 | 0.007 | 0.289 | ||
| fixed | flock size (ln) | -0.111 | -0.217 | -0.006 | ||
| fixed | body mass (ln) | 0.32 | 0.164 | 0.481 | ||
| fixed | time (sine of radians) | 0.117 | -0.027 | 0.261 | ||
| fixed | time (cosine of radians) | 0.031 | -0.099 | 0.15 | ||
| fixed | temperaturre | -0.018 | -0.155 | 0.115 | ||
| fixed | Google Mobility | -0.088 | -0.232 | 0.052 | ||
| random | species within site (Intercept) | 15% | 14% | 15% | ||
| random | species within day & year (Intercept) | 3% | 3% | 3% | ||
| random | site (Intercept) | 17% | 14% | 20% | ||
| random | species (Intercept) | 0% | 0% | 0% | ||
| random | genus (Intercept) | 3% | 1% | 4% | ||
| random | weekday (Intercept) | 0% | 0% | 0% | ||
| random | Residual | 62% | 57% | 67% | ||
| Poland | 329 | fixed | (Intercept) | 0.028 | -0.198 | 0.259 |
| fixed | year | -0.25 | -0.407 | -0.088 | ||
| fixed | starting distance (ln)) | 0.532 | 0.447 | 0.617 | ||
| fixed | flock size (ln) | -0.03 | -0.112 | 0.05 | ||
| fixed | body mass (ln) | -0.082 | -0.19 | 0.024 | ||
| fixed | time (sine of radians) | 0.135 | -0.135 | 0.402 | ||
| fixed | time (cosine of radians) | -0.252 | -0.517 | 0.013 | ||
| fixed | temperaturre | -0.201 | -0.306 | -0.1 | ||
| fixed | Google Mobility | 0 | -0.14 | 0.14 | ||
| random | species within day & year (Intercept) | 7% | 6% | 6% | ||
| random | species (Intercept) | 9% | 7% | 10% | ||
| random | genus (Intercept) | 0% | 0% | 0% | ||
| random | weekday (Intercept) | 11% | 5% | 16% | ||
| random | Residual | 74% | 67% | 82% | ||
| Poland, without starting distance | 329 | fixed | (Intercept) | 0.018 | -0.242 | 0.271 |
| fixed | year | -0.4 | -0.584 | -0.225 | ||
| fixed | flock size (ln) | 0.052 | -0.042 | 0.149 | ||
| fixed | body mass (ln) | -0.008 | -0.159 | 0.137 | ||
| fixed | time (sine of radians) | -0.02 | -0.328 | 0.284 | ||
| fixed | time (cosine of radians) | -0.087 | -0.406 | 0.215 | ||
| fixed | temperaturre | -0.189 | -0.311 | -0.071 | ||
| fixed | Google Mobility | 0.043 | -0.114 | 0.203 | ||
| random | species within day & year (Intercept) | 5% | 4% | 5% | ||
| random | species (Intercept) | 21% | 17% | 25% | ||
| random | genus (Intercept) | 0% | 0% | 0% | ||
| random | weekday (Intercept) | 5% | 2% | 8% | ||
| random | Residual | 70% | 63% | 77% | ||
| Czechia | 1054 | fixed | (Intercept) | 0.229 | 0.025 | 0.438 |
| fixed | starting distance (ln)) | 0.344 | 0.286 | 0.403 | ||
| fixed | flock size (ln) | 0.01 | -0.04 | 0.061 | ||
| fixed | body mass (ln) | 0.184 | 0.025 | 0.348 | ||
| fixed | time (sine of radians) | 0.063 | -0.008 | 0.133 | ||
| fixed | time (cosine of radians) | 0.034 | -0.022 | 0.09 | ||
| fixed | temperaturre | 0.066 | -0.005 | 0.137 | ||
| fixed | Google Mobility | -0.09 | -0.162 | -0.018 | ||
| random | species within site (Intercept) | 5% | 5% | 5% | ||
| random | species within day & year (Intercept) | 1% | 0% | 1% | ||
| random | species (Intercept) | 14% | 11% | 16% | ||
| random | genus (Intercept) | 14% | 10% | 18% | ||
| random | site (Intercept) | 6% | 4% | 7% | ||
| random | weekday (Intercept) | 0% | 0% | 0% | ||
| random | Residual | 61% | 55% | 69% | ||
| Czechia, without starting distance | 1054 | fixed | (Intercept) | 0.227 | 0.02 | 0.438 |
| fixed | starting distance (ln)) | 0.344 | 0.281 | 0.402 | ||
| fixed | flock size (ln) | 0.01 | -0.039 | 0.059 | ||
| fixed | body mass (ln) | 0.184 | 0.031 | 0.339 | ||
| fixed | time (sine of radians) | 0.063 | -0.008 | 0.13 | ||
| fixed | time (cosine of radians) | 0.036 | -0.021 | 0.092 | ||
| fixed | temperaturre | 0.065 | -0.007 | 0.134 | ||
| fixed | Google Mobility | -0.09 | -0.161 | -0.017 | ||
| random | species within site (Intercept) | 5% | 5% | 5% | ||
| random | species within day & year (Intercept) | 1% | 0% | 1% | ||
| random | species (Intercept) | 14% | 11% | 16% | ||
| random | genus (Intercept) | 14% | 10% | 18% | ||
| random | site (Intercept) | 6% | 4% | 7% | ||
| random | weekday (Intercept) | 0% | 0% | 0% | ||
| random | Residual | 61% | 55% | 69% | ||
| Hungary | 874 | fixed | (Intercept) | 0.074 | -0.211 | 0.371 |
| fixed | year | 0.167 | 0.055 | 0.278 | ||
| fixed | starting distance (ln)) | 0.502 | 0.442 | 0.564 | ||
| fixed | flock size (ln) | -0.001 | -0.053 | 0.051 | ||
| fixed | body mass (ln) | 0.001 | -0.142 | 0.141 | ||
| fixed | time (sine of radians) | -0.126 | -0.193 | -0.059 | ||
| fixed | time (cosine of radians) | 0.017 | -0.052 | 0.083 | ||
| fixed | temperaturre | -0.033 | -0.111 | 0.045 | ||
| fixed | Google Mobility | -0.027 | -0.133 | 0.077 | ||
| random | species within day & year (Intercept) | 2% | 2% | 2% | ||
| random | species within site (Intercept) | 0% | 0% | 0% | ||
| random | species (Intercept) | 1% | 1% | 1% | ||
| random | genus (Intercept) | 12% | 8% | 16% | ||
| random | weekday (Intercept) | 1% | 0% | 2% | ||
| random | site (Intercept) | 9% | 5% | 13% | ||
| random | Residual | 74% | 65% | 83% | ||
| Hungary, without starting distance | 874 | fixed | (Intercept) | 0.062 | -0.411 | 0.527 |
| fixed | year | 0.314 | 0.188 | 0.441 | ||
| fixed | flock size (ln) | -0.005 | -0.064 | 0.056 | ||
| fixed | body mass (ln) | 0.287 | 0.066 | 0.503 | ||
| fixed | time (sine of radians) | -0.166 | -0.242 | -0.088 | ||
| fixed | time (cosine of radians) | 0.035 | -0.035 | 0.107 | ||
| fixed | temperaturre | -0.049 | -0.145 | 0.045 | ||
| fixed | Google Mobility | -0.146 | -0.269 | -0.022 | ||
| random | species within day & year (Intercept) | 2% | 2% | 2% | ||
| random | species within site (Intercept) | 1% | 1% | 1% | ||
| random | species (Intercept) | 1% | 1% | 2% | ||
| random | genus (Intercept) | 26% | 20% | 32% | ||
| random | weekday (Intercept) | 2% | 1% | 3% | ||
| random | site (Intercept) | 16% | 13% | 19% | ||
| random | Residual | 52% | 42% | 63% | ||
| Australia | 1065 | fixed | (Intercept) | -0.053 | -0.192 | 0.097 |
| fixed | year | 0.028 | -0.035 | 0.091 | ||
| fixed | starting distance (ln)) | 0.5 | 0.446 | 0.552 | ||
| fixed | flock size (ln) | 0.025 | -0.022 | 0.072 | ||
| fixed | body mass (ln) | -0.049 | -0.148 | 0.045 | ||
| fixed | time (sine of radians) | -0.044 | -0.117 | 0.03 | ||
| fixed | time (cosine of radians) | -0.044 | -0.111 | 0.026 | ||
| fixed | temperaturre | 0.017 | -0.039 | 0.076 | ||
| fixed | Google Mobility | -0.036 | -0.093 | 0.02 | ||
| random | species within day & year (Intercept) | 6% | 6% | 6% | ||
| random | species within site (Intercept) | 12% | 12% | 12% | ||
| random | species (Intercept) | 11% | 9% | 13% | ||
| random | genus (Intercept) | 2% | 2% | 2% | ||
| random | site (Intercept) | 8% | 6% | 10% | ||
| random | weekday (Intercept) | 0% | 0% | 0% | ||
| random | Residual | 61% | 56% | 65% | ||
| Australia, without starting distance | 1065 | fixed | (Intercept) | -0.103 | -0.296 | 0.095 |
| fixed | year | 0.087 | 0.013 | 0.163 | ||
| fixed | flock size (ln) | 0.054 | -0.001 | 0.11 | ||
| fixed | body mass (ln) | 0.066 | -0.058 | 0.188 | ||
| fixed | time (sine of radians) | -0.028 | -0.114 | 0.057 | ||
| fixed | time (cosine of radians) | -0.019 | -0.099 | 0.062 | ||
| fixed | temperaturre | 0.01 | -0.059 | 0.079 | ||
| fixed | Google Mobility | -0.027 | -0.097 | 0.041 | ||
| random | species within day & year (Intercept) | 7% | 7% | 7% | ||
| random | species within site (Intercept) | 11% | 11% | 11% | ||
| random | species (Intercept) | 13% | 11% | 15% | ||
| random | genus (Intercept) | 6% | 5% | 8% | ||
| random | site (Intercept) | 10% | 8% | 12% | ||
| random | weekday (Intercept) | 0% | 0% | 0% | ||
| random | Residual | 52% | 48% | 57% |
Species trends for stringency or Google Mobility
ss <- s[Nsp > 9]
ss[, sp2 := gsub(" ", "\n", sp)]
ss[, Country := factor(Country, levels = rev(c("Finland", "Poland", "Czechia", "Hungary", "Australia")))]
gs2 <-
ggplot(ss, aes(x = StringencyIndex, y = FID)) +
# stat_smooth(method = 'rlm', se = FALSE, col = 'black', lwd = 0.5)+
geom_point(pch = 21, alpha = 0.7, aes(fill = Country), col = 'grey20', alpha =0.8)+
stat_smooth(se = FALSE, aes(colour = "Locally weighted\nsmoothing"), lwd = 0.5) + # show_guide=TRUE
facet_wrap(~sp2, ncol = 6) +
scale_fill_manual(values = col3__,guide = guide_legend(reverse = TRUE)) +
#scale_fill_viridis(discrete = TRUE, guide = guide_legend(reverse = FALSE)) +
scale_x_continuous("Stringency index of governmental COVID-19 restrictions", expand = c(0, 0)) +
scale_y_continuous("Flight initiation distance [m]", expand = c(0, 0), trans = "log10") +
# annotate("text", x = 1, y = 1, label = c(rep("", 52),"Observation"), hjust = -0.08, size = 1) +
# labs(title = "Species means per sampling location")+
scale_colour_manual(values = c("grey60")) +
# scale_color_manual(name = 'try', values = c('LOESS smoothed = "grey60"'))+
theme_MB +
theme(
plot.title = element_text(size = 7),
strip.background = element_blank(),
# strip.text.x = element_text(size = 4.5, color="grey30", margin=margin(1,1,1,1,"mm")),
# panel.spacing = unit(1, "mm"),
legend.position = c(1, 0.01),
legend.justification = c(1, 0),
legend.title = element_blank(),
# legend.spacing.y = unit(-0.78, "cm")
# legend.spacing.y = unit(0.02, "cm") use if LOESS smooth text as legend
legend.spacing.y = unit(-0.9, "cm")
)
gsg2 <- ggplotGrob(gs2) # gg$layout$name
gsgx2 <- gtable_filter_remove(gsg2, name = paste0("axis-b-", c(2, 4), "-9"), trim = FALSE)
grid.draw(gsgx2)
if(save_plot==TRUE){
ggsave(here::here("Outputs/Fig_S7_width-152mm.png"), gsgx2, width = 15.24, height = 19, unit = "cm", dpi = 600)
}Figure S7 | Species-specific avian tolerance towards humans in relation to stringency of antipandemic governmental restrictions during COVID-19 shutdowns quantified as a stringency index. Each dot represents a single escape distance observation (not corrected for other factors such as starting distance of the observer) and a day-specific value of governmental stringency index in a given country. Dot colour highlights the country. Grey lines represent locally weighted smoothing, a non-parametric local regression fitted with the ggplot function of ggplot2 package (Wickham 2016), highlighting heterogenous (and usually unclear – close to zero) within- and between- species trends. Note, the y-axes is on the log-scale, some species lack trend lines because data distribution hindered the smoothing and visualised are only data for species with ≥10 escape distance observations.
ss[, NspC := .N, by = "sp_country"]
ssc <- ss[NspC > 9]
ssc[, sp2 := gsub(" ", "\n", sp)]
ssc[, Google_mobility := parks_percent_change_from_baseline]
ssc[, Country := factor(Country, levels = rev(c("Finland", "Poland", "Czechia", "Hungary", "Australia")))]
# two rows labels
gt2 <-
ggplot(ssc, aes(x = Google_mobility, y = FID, )) +
# stat_smooth(method = 'rlm', se = FALSE, col = 'black', lwd = 0.5)+
geom_point(pch = 21, alpha = 0.7, aes(fill = Country), col = 'grey20', alpha =0.8)+
stat_smooth(se = FALSE, aes(colour = "Locally weighted\nsmoothing"), lwd = 0.5) + # show_guide=TRUE
facet_wrap(~sp2, ncol = 6) +
scale_fill_manual(values = col3__,guide = guide_legend(reverse = TRUE)) +
scale_colour_manual(values = c("grey60")) +
scale_x_continuous("Google Mobility", expand = c(0, 0)) +
scale_y_continuous("Flight initiation distance [m]", expand = c(0, 0), trans = "log10") +
# annotate("text", x = 1, y = 1, label = c(rep("", 52),"Observation"), hjust = -0.08, size = 1) +
# labs(title = "Species means per sampling location")+
# scale_colour_manual(values=c('grey60'))+
# scale_color_manual(name = 'try', values = c('LOESS smoothed = "grey60"'))+
theme_MB +
theme(
plot.title = element_text(size = 7),
strip.background = element_blank(),
# strip.text.x = element_text(size = 4.5, color="grey30", margin=margin(1,1,1,1,"mm")),
# panel.spacing = unit(1, "mm"),
legend.position = c(1, 0.01),
legend.justification = c(1, 0),
legend.title = element_blank(),
# legend.spacing.y = unit(-0.78, "cm")
# legend.spacing.y = unit(0.02, "cm") use if LOESS smooth text as legend
legend.spacing.y = unit(-0.9, "cm")
)
gtg2 <- ggplotGrob(gt2) # gg$layout$name
gtgx2 <- gtable_filter_remove(gtg2, name = c("axis-b-2-9", "axis-b-5-8"), trim = FALSE) # paste0("axis-b-", c(2, 4), "-9")
grid.draw(gtgx2)
if(save_plot==TRUE){
ggsave(here::here("Outputs/Fig_S7_width-152mm_Google.png"), gtgx2, width = 15.24, height = 19, unit = "cm", dpi = 600)
}Figure S8 | Species-specific avian tolerance towards humans in relation to Google Mobility. Each dot represents a single escape distance observation (not corrected for other factors such as starting distance of the observer) and a day-specific value of Google Mobility for parks in a given country. Dot colour highlights the country. Grey lines represent locally weighted smoothing, a non-parametric local regression fitted with the ggplot function of ggplot2 package (Wickham 2016), highlighting heterogenous (and usually unclear – close to zero) within- and between- species trends. The y-axes is on the log-scale. Some species lack trend lines because data distribution hindered the smoothing and visualised are only data for species with ≥10 escape distance observations.
Testing for phylo-signal in residuals
Example code. Don’t run, as it takes too long. To make the Table S5, load the saved mcmc outputs instead.
# packages and data
library(ape)
library(brms)
library(coda)
library(dplyr)
library(geiger)
library(parallel)
library(phangorn)
library(phylobase)
sapply(c(
"magrittr", "car",
"ape", "phytools", "stringr", "plyr"
),
require,
character.only = TRUE
)
trees = read.tree("Data/trees.tre")
d[, scinam := Species]
s[, scinam := Species]
# get residuals
# period
ms <- lmer(scale(log(FID)) ~
scale(log(SD)) +
scale(log(FlockSize)) +
scale(log(BodyMass)) +
scale(sin(rad)) + scale(cos(rad)) +
# scale(Day)+
scale(Temp) +
scale(Covid) +
(1 | Year) + (1 | weekday) + (1 | genus) + (1 | Species) + (1 | sp_day_year) + (scale(Covid) | Country) + (scale(Covid) | IDLocality) + (1 | sp_loc),
data = d, REML = FALSE,
control <- lmerControl(
optimizer = "optimx", optCtrl = list(method = "nlminb")
)
)
d[, res_period := resid(ms)]
# stringency
m01a <- lmer(scale(log(FID)) ~
scale(Year) +
scale(log(SD)) +
scale(log(FlockSize)) +
scale(log(BodyMass)) +
scale(sin(rad)) + scale(cos(rad)) +
# scale(Day)+
scale(Temp) +
scale(StringencyIndex) +
(1 | weekday) + (1 | genus) + (1 | Species) + (1 | sp_day_year) + (scale(StringencyIndex) | Country) + (1 | IDLocality) + (1 | sp_loc),
data = s, REML = FALSE
)
s[, res := resid(m01a)]
# google
ss = s[!is.na(parks_percent_change_from_baseline)]
ss[, country_year := paste(Country, Year)] # table(paste(s$Country, s$Year))
ss[parks_percent_change_from_baseline < 0, google := "before_zero"]
ss[parks_percent_change_from_baseline > 0, google := "after_zero"]
ss[, sp_country_google := paste(sp_country, google)]
mg01a = lmer(scale(log(FID)) ~
scale(Year) +
scale(log(SD)) +
scale(log(FlockSize)) +
scale(log(BodyMass)) +
scale(sin(rad)) + scale(cos(rad)) +
# scale(Day)+
scale(Temp) +
scale(parks_percent_change_from_baseline) +
(1 | weekday) + (1 | genus) + (1 | Species) + (1 | sp_day_year) + (scale(parks_percent_change_from_baseline) | Country) + (1 | IDLocality) + (1 | sp_loc),
data = ss, REML = FALSE
)
ss[, res := resid(mg01a)]
# MCC tree and covariance matrix
namesd <- d %>% distinct(Species)
row.names(namesd) <- namesd$Species
namecheck_fidd <- name.check(trees[[1]], namesd)
trees_fidd <- lapply(trees, drop.tip, namecheck_fidd$tree_not_data)
class(trees_fidd) <- "multiPhylo"
length(trees_fidd[[1]]$tip.label)
tree_fidd <- maxCladeCred(trees_fidd)
tree_fidd$node.label <- NULL
inv.phylo_d <- inverseA(tree_fidd, nodes = "ALL", scale = TRUE)
phyloresd <- solve(inv.phylo_d$Ainv)
rownames(phyloresd) <- rownames(inv.phylo_d$Ainv)
namess <- s %>% distinct(Species)
row.names(namess) <- namess$Species
namecheck_fids <- name.check(trees[[1]], namess)
trees_fids <- lapply(trees, drop.tip, namecheck_fids$tree_not_data)
class(trees_fids) <- "multiPhylo"
length(trees_fids[[1]]$tip.label)
tree_fids <- maxCladeCred(trees_fids)
tree_fids$node.label <- NULL
inv.phylo_s <- inverseA(tree_fids, nodes = "ALL", scale = TRUE)
phyloress <- solve(inv.phylo_s$Ainv)
rownames(phyloress) <- rownames(inv.phylo_s$Ainv)
namesss <- ss %>% distinct(Species)
row.names(namesss) <- namesss$Species
namecheck_fidss <- name.check(trees[[1]], namesss)
trees_fidss <- lapply(trees, drop.tip, namecheck_fidss$tree_not_data)
class(trees_fidss) <- "multiPhylo"
length(trees_fidss[[1]]$tip.label)
tree_fidss <- maxCladeCred(trees_fidss)
tree_fidss$node.label <- NULL
inv.phylo_ss <- inverseA(tree_fidss, nodes = "ALL", scale = TRUE)
phyloresss <- solve(inv.phylo_ss$Ainv)
rownames(phyloresss) <- rownames(inv.phylo_ss$Ainv)
# brms
# period - no - need to run (takes time), if you load the above DAT_brms.Rdata
# without phylo
priors <- get_prior(res_period ~ 0 + Intercept + (1 | Species), data = d) #
mP_no = brm(
form = res_period ~ 0 + Intercept + (1 | Species),
data = d,
cores = 2,
chains = 5,
control = list(adapt_delta = 0.99),
iter = 5000,
thin = 5,
sample_prior = "yes",
save_pars = save_pars(all = TRUE),
prior = priors,
seed = 5
)
plot(mP_no, ask = FALSE)
pp_check(mP_no, ndraws = 100)
mcmc_plot(mP_no, type = "acf")
summary(mP_no)
# with phylo
priors <- get_prior(res_period ~ 0 + Intercept + (1 | Species) + (1 | gr(scinam, cov = A)), data = d, data2 = list(A = phyloresd))
mP_yes = brm(
form = res_period ~ 0 + Intercept + (1 | Species) + (1 | gr(scinam, cov = A)),
data = d,
data2 = list(A = phyloresd),
cores = 2,
chains = 5,
control = list(adapt_delta = 0.99),
iter = 5000,
thin = 5,
sample_prior = "yes",
save_pars = save_pars(all = TRUE),
prior = priors,
seed = 5
)
plot(mP_yes, ask = FALSE)
pp_check(mP_yes, ndraws = 100)
# pp_check(m_yes, type = "scatter_avg_grouped", group = "Species") + geom_abline(intercept = 0, slope = 1 , color = "red", lty = 2)
mcmc_plot(mP_yes, type = "acf")
summary(mP_yes)
# lambda
hypothesis(mP_yes, "sd_scinam__Intercept^2 / (sd_scinam__Intercept^2 + sd_Species__Intercept^2) = 0", class = NULL)
hypothesis(mP_yes, "sd_scinam__Intercept^2 / (sd_scinam__Intercept^2 + sd_Species__Intercept^2 + sigma^2) = 0", class = NULL)
v_sc <- (VarCorr(mP_yes, summary = FALSE)$scinam$sd)^2
v_sp <- (VarCorr(mP_yes, summary = FALSE)$Species$sd)^2
v_r <- (VarCorr(mP_yes, summary = FALSE)$residual$sd)^2
summary(as.mcmc(v_sc / (v_sc + v_sp + v_r)))
# compare period models
# 1. LOOic
loo(mP_yes)
loo(mP_no)
# 2. Bayes factor
bayes_factor(mP_no, mP_yes)
# 3. Poterior probability
post_prob(mP_yes, mP_no)
# stringency index - no - need to run (takes time), if you load the above DAT_brms.Rdata
# without phylo
priors <- get_prior(res ~ 0 + Intercept + (1 | Species), data = s) #
m_no = brm(
form = res ~ 0 + Intercept + (1 | Species),
data = s,
cores = 2,
chains = 5,
control = list(adapt_delta = 0.99),
iter = 5000,
thin = 5,
sample_prior = "yes",
save_pars = save_pars(all = TRUE),
prior = priors,
seed = 5
)
plot(m_no, ask = FALSE)
pp_check(m_no, ndraws = 100)
mcmc_plot(m_no, type = "acf")
summary(m_no)
# with phylo
priors <- get_prior(res ~ 0 + Intercept + (1 | Species) + (1 | gr(scinam, cov = A)), data = s, data2 = list(A = phyloress))
m_yes = brm(
form = res ~ 0 + Intercept + (1 | Species) + (1 | gr(scinam, cov = A)),
data = s,
data2 = list(A = phyloress),
cores = 2,
chains = 5,
control = list(adapt_delta = 0.99),
iter = 5000,
thin = 5,
sample_prior = "yes",
save_pars = save_pars(all = TRUE),
prior = priors,
seed = 5
)
plot(m_yes, ask = FALSE)
pp_check(m_yes, ndraws = 100)
# pp_check(m_yes, type = "scatter_avg_grouped", group = "Species") + geom_abline(intercept = 0, slope = 1 , color = "red", lty = 2)
mcmc_plot(m_yes, type = "acf")
summary(m_yes)
# lambda
hypothesis(m_yes, "sd_scinam__Intercept^2 / (sd_scinam__Intercept^2 + sd_Species__Intercept^2) = 0", class = NULL) # lambda
hypothesis(m_yes, "sd_scinam__Intercept^2 / (sd_scinam__Intercept^2 + sd_Species__Intercept^2 + sigma^2) = 0", class = NULL)
v_sc <- (VarCorr(m_yes, summary = FALSE)$scinam$sd)^2
v_sp <- (VarCorr(m_yes, summary = FALSE)$Species$sd)^2
v_r <- (VarCorr(m_yes, summary = FALSE)$residual$sd)^2
summary(as.mcmc(v_sc / (v_sc + v_sp + v_r)))
Mode <- function(x) {
ux <- unique(x)
ux[which.max(tabulate(match(x, ux)))]
}
Mode(as.mcmc(v_sc / (v_sc + v_sp + v_r)))
# compare stringency models
# 1. LOOic
loo(m_yes)
loo(m_no)
# 2. Bayes factor
bayes_factor(m_no, m_yes)
# 3. Poterior probability
post_prob(m_yes, m_no)
# Google Mobility - no - need to run (takes time), if you load the above DAT_brms.Rdata
# without phylo
priors <- get_prior(res ~ 0 + Intercept + (1 | Species), data = ss) #
mG_no = brm(
form = res ~ 0 + Intercept + (1 | Species),
data = ss,
cores = 2,
chains = 5,
control = list(adapt_delta = 0.99),
iter = 5000,
thin = 5,
sample_prior = "yes",
save_pars = save_pars(all = TRUE),
prior = priors,
seed = 5
)
plot(mG_no, ask = FALSE)
pp_check(mG_no, ndraws = 100)
mcmc_plot(mG_no, type = "acf")
summary(mG_no)
# with phylo
priors <- get_prior(res ~ 0 + Intercept + (1 | Species) + (1 | gr(scinam, cov = A)), data = ss, data2 = list(A = phyloress))
mG_yes = brm(
form = res ~ 0 + Intercept + (1 | Species) + (1 | gr(scinam, cov = A)),
data = ss,
data2 = list(A = phyloress),
cores = 2,
chains = 5,
control = list(adapt_delta = 0.99),
iter = 5000,
thin = 5,
sample_prior = "yes",
save_pars = save_pars(all = TRUE),
prior = priors,
seed = 5
)
plot(mG_yes, ask = FALSE)
pp_check(mG_yes, ndraws = 100)
# pp_check(m_yes, type = "scatter_avg_grouped", group = "Species") + geom_abline(intercept = 0, slope = 1 , color = "red", lty = 2)
mcmc_plot(mG_yes, type = "acf")
summary(mG_yes)
# lambda
hypothesis(mG_yes, "sd_scinam__Intercept^2 / (sd_scinam__Intercept^2 + sd_Species__Intercept^2) = 0", class = NULL) # lambda
hypothesis(mG_yes, "sd_scinam__Intercept^2 / (sd_scinam__Intercept^2 + sd_Species__Intercept^2 + sigma^2) = 0", class = NULL)
v_sc <- (VarCorr(mG_yes, summary = FALSE)$scinam$sd)^2
v_sp <- (VarCorr(mG_yes, summary = FALSE)$Species$sd)^2
v_r <- (VarCorr(mG_yes, summary = FALSE)$residual$sd)^2
summary(as.mcmc(v_sc / (v_sc + v_sp + v_r)))
Mode <- function(x) {
ux <- unique(x)
ux[which.max(tabulate(match(x, ux)))]
}
Mode(as.mcmc(v_sc / (v_sc + v_sp + v_r)))
# compare stringency models
# 1. LOOic
loo(mG_yes)
loo(mG_no)
# 2. Bayes factor
bayes_factor(mG_no, mG_yes)
# 3. Poterior probability
post_prob(mG_yes, mG_no)
# export models
save(file = 'Data/rev_DAT_brms.Rdata', mP_no, mP_yes, m_no, m_yes, mG_no, mG_yes)
load(here::here("Data/rev_DAT_brms.Rdata"))
v_sc <- (VarCorr(mP_yes, summary = FALSE)$scinam$sd)^2
v_sp <- (VarCorr(mP_yes, summary = FALSE)$Species$sd)^2
v_r <- (VarCorr(mP_yes, summary = FALSE)$residual$sd)^2
s_v_sc <- (VarCorr(m_yes, summary = FALSE)$scinam$sd)^2
s_v_sp <- (VarCorr(m_yes, summary = FALSE)$Species$sd)^2
s_v_r <- (VarCorr(m_yes, summary = FALSE)$residual$sd)^2
g_v_sc <- (VarCorr(mG_yes, summary = FALSE)$scinam$sd)^2
g_v_sp <- (VarCorr(mG_yes, summary = FALSE)$Species$sd)^2
g_v_r <- (VarCorr(mG_yes, summary = FALSE)$residual$sd)^2
# summary(as.mcmc(v_sc / (v_sc + v_sp + v_r)))
ts5 <- data.table(
"Residuals\nfrom" = c(
"Table S2a - model 01a",
"Table S2b - model 01a",
"Table S2c - model 01a"
),
"Variance explained by phylogeny\n(95%CI)" = c(
paste0(paste(round(quantile(as.mcmc(v_sc / (v_sc + v_sp + v_r)), probs = c(0.025, 0.975)) * 100, 2), collapse = "-"), "%"),
paste0(paste(round(quantile(as.mcmc(s_v_sc / (s_v_sc + s_v_sp + s_v_r)), probs = c(0.025, 0.975)) * 100, 2), collapse = "-"), "%"),
paste0(paste(round(quantile(as.mcmc(g_v_sc / (g_v_sc + g_v_sp + g_v_r)), probs = c(0.025, 0.975)) * 100, 2), collapse = "-"), "%")
),
"Bayes\nfactor" = c(
round(bayes_factor(mP_no, mP_yes)$bf),
round(bayes_factor(m_no, m_yes)$bf),
round(bayes_factor(mG_no, mG_yes)$bf)
),
"Probability of non-phylogenetic model" = c(
round(post_prob(mP_yes, mP_no)[2], 2),
round(post_prob(m_yes, m_no)[2], 2),
round(post_prob(mG_yes, mG_no)[2], 2)
)
)
save(file = "Data/T_S5.Rdata", ts5)
fwrite(file = "Data/T_S5.csv", ts5)Table S5 | Comparison of models without and with control for phylogeny.
ts5 = fread(here::here('Data/T_S5.csv'))
ts5 %>%
kbl(align=c('l', 'r', 'r','r')) %>%
kable_paper("hover", full_width = F)| Residuals from | Variance explained by phylogeny (95%CI) | Bayes factor | Probability of non-phylogenetic model |
|---|---|---|---|
| Table S2a - model 01a | 0-0.19% | 99 | 0.98 |
| Table S2b - model 01a | 0-0.38% | 115 | 0.99 |
| Table S2c - model 01a | 0-0.36% | 96 | 0.99 |
Variance explained by phylogeny (95%CI) represents percentage of variance explained by phylogeny in an intercept only model fitted to residuals of the original models in STAN (Stan Development Team 2022) using ‘brm’ function from ‘brms’ R-package (Bürkner 2017, Bürkner 2018) with phylogeny and species as random effects. Bayes factor in favour of model without phylogeny and probability of non-phylogenetic model in comparison to a model with phylogeny. The Gelman-Rubin diagnostics was 1 for all models, indicating model convergence (Brooks & Gelman 1998). Note that for all cases, the model without phylogeny fits residuals better than a model with phylogeny, which justifies our use of simple original models, not controlled for phylogeny.
# ENDReferences
- Bulla, M., Blumstein, D.T., Benedetti, Y., Floigl, K., Jokimäki, J., Kaisanlahti-Jokimäki, M.-L., Markó, G., Morelli, F., Siretckaia, A., Szakony, S., Weston, M.A., Zeid, F.A., Tryjanowski, P., Albrecht, T. & Mikula, P. (2022). Supporting information for ‘Urban birds’ flight responses were unaffected by the COVID-19 shutdowns’. Open Science Framework https://doi.org/10.17605/OSF.IO/WUZH7.
- Bürkner PC. (2018). Advanced Bayesian multilevel modeling with the R package brms. R Journal 10:395–411. https://doi.org/10.32614/RJ-2018-017
- Bürkner PC. (2017). brms: An R package for Bayesian multilevel models using Stan. J Stat Softw 80:1–28. https://doi.org/10.18637/jss.v080.i01
- Gelman, A., Su, Y.-S., Yajima, M., Hill, J., Pittau, M., Kerman, J., Zheng, T., & Vincent, D. (2016). Data Analysis using Regression and Multilevel/Hierarchical Models. In CRAN Repository (1.8-6.; pp. 1–53).
- Hale T., Angrist N., Goldszmidt R., Kira B., Petherick A., Phillips T., Webster S., Cameron-Blake E., Hallas L., Majumdar S., Tatlow H. (2021). A global panel database of pandemic policies (Oxford COVID-19 Government Response Tracker). Nature Human Behaviour 2021 5:4 5:529–538. https://doi.org/10.1038/s41562-021-01079-8.
- Peterson, B. G., & Carl, P. (2020). PerformanceAnalytics: Econometric Tools for Performance and Risk Analysis. R package version 2.0.4. https://CRAN.R-project.org/package=PerformanceAnalytics
- Stan Development Team (2022). Stan Modeling Language Users Guide and Reference Manual, Version 2.28. https://mc-stan.org/users/documentation/
- Wickham, H. (2016). ggplot2: Elegant Graphics for Data Analysis. Springer-Verlag, New York.
Session info
sessionInfo()## R version 4.2.0 (2022-04-22)
## Platform: x86_64-apple-darwin17.0 (64-bit)
## Running under: macOS Mojave 10.14.6
##
## Matrix products: default
## BLAS: /Library/Frameworks/R.framework/Versions/4.2/Resources/lib/libRblas.0.dylib
## LAPACK: /Library/Frameworks/R.framework/Versions/4.2/Resources/lib/libRlapack.dylib
##
## locale:
## [1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
##
## attached base packages:
## [1] parallel grid stats graphics grDevices utils datasets methods base
##
## other attached packages:
## [1] rmarkdown_2.18 plyr_1.8.7 stringr_1.5.0 phytools_1.3-1 maps_3.4.1 car_3.0-13
## [7] magrittr_2.0.3 phylobase_0.8.10 phangorn_2.10.0 MCMCglmm_2.33 geiger_2.0.9 dplyr_1.0.10
## [13] coda_0.19-4 brms_2.17.0 Rcpp_1.0.9 ape_5.6-2 viridis_0.6.2 viridisLite_0.4.1
## [19] scales_1.2.1 rphylopic_0.3.0 rmeta_3.0 RColorBrewer_1.1-3 png_0.1-7 PerformanceAnalytics_2.0.4
## [25] xts_0.12.1 zoo_1.8-10 performance_0.9.0 optimx_2022-4.30 multcomp_1.4-19 TH.data_1.1-1
## [31] survival_3.3-1 mvtnorm_1.1-3 kableExtra_1.3.4 here_1.0.1 gtable_0.3.1 ggtext_0.1.2
## [37] ggsci_2.9 ggpubr_0.4.0 ggimage_0.3.1 ggplot2_3.3.6 foreach_1.5.2 effects_4.2-1
## [43] carData_3.0-5 data.table_1.14.6 arm_1.12-2 lme4_1.1-29 Matrix_1.4-1 MASS_7.3-57
##
## loaded via a namespace (and not attached):
## [1] utf8_1.2.2 tidyselect_1.1.2 htmlwidgets_1.5.4 combinat_0.0-8 RNeXML_2.4.7 munsell_0.5.0 codetools_0.2-18
## [8] DT_0.23 miniUI_0.1.1.1 withr_2.5.0 Brobdingnag_1.2-7 colorspace_2.0-3 highr_0.9 uuid_1.1-0
## [15] knitr_1.41 rstudioapi_0.14 stats4_4.2.0 ggsignif_0.6.3 bayesplot_1.9.0 labeling_0.4.2 emmeans_1.7.4-1
## [22] rstan_2.21.5 optimParallel_1.0-2 mnormt_2.1.1 farver_2.1.1 bridgesampling_1.1-2 rprojroot_2.0.3 vctrs_0.5.1
## [29] generics_0.1.3 clusterGeneration_1.3.7 xfun_0.35 R6_2.5.1 markdown_1.1 cachem_1.0.6 gridGraphics_0.5-1
## [36] assertthat_0.2.1 promises_1.2.0.1 nnet_7.3-17 processx_3.8.0 sandwich_3.0-1 rlang_1.0.6 systemfonts_1.0.4
## [43] scatterplot3d_0.3-42 splines_4.2.0 lazyeval_0.2.2 rstatix_0.7.0 broom_0.8.0 checkmate_2.1.0 inline_0.3.19
## [50] yaml_2.3.6 reshape2_1.4.4 abind_1.4-5 threejs_0.3.3 crosstalk_1.2.0 backports_1.4.1 httpuv_1.6.5
## [57] gridtext_0.1.4 tensorA_0.36.2 tools_4.2.0 ggplotify_0.1.0 gridBase_0.4-7 cubature_2.0.4.4 ellipsis_0.3.2
## [64] jquerylib_0.1.4 posterior_1.2.2 ggridges_0.5.3 progress_1.2.2 base64enc_0.1-3 purrr_0.3.4 ps_1.7.2
## [71] prettyunits_1.1.1 deSolve_1.32 survey_4.1-1 crul_1.2.0 magick_2.7.3 colourpicker_1.1.1 matrixStats_0.62.0
## [78] hms_1.1.2 shinyjs_2.1.0 mime_0.12 evaluate_0.18 xtable_1.8-4 XML_3.99-0.9 shinystan_2.6.0
## [85] gridExtra_2.3 rstantools_2.2.0 compiler_4.2.0 tibble_3.1.8 crayon_1.5.1 minqa_1.2.4 StanHeaders_2.21.0-7
## [92] htmltools_0.5.4 corpcor_1.6.10 ggfun_0.0.6 later_1.3.0 tidyr_1.2.0 expm_0.999-6 RcppParallel_5.1.5
## [99] DBI_1.1.3 subplex_1.8 boot_1.3-28 ade4_1.7-19 cli_3.4.1 mitools_2.4 quadprog_1.5-8
## [106] insight_0.18.8 igraph_1.3.5 pkgconfig_2.0.3 rncl_0.8.6 numDeriv_2016.8-1.1 xml2_1.3.3 dygraphs_1.1.1.6
## [113] svglite_2.1.0 bslib_0.4.1 webshot_0.5.4 estimability_1.3 rvest_1.0.3 yulab.utils_0.0.4 distributional_0.3.0
## [120] callr_3.7.3 digest_0.6.30 httpcode_0.3.0 fastmatch_1.1-3 curl_4.3.3 shiny_1.7.1 gtools_3.9.2
## [127] nloptr_2.0.3 lifecycle_1.0.3 nlme_3.1-157 jsonlite_1.8.4 fansi_1.0.3 pillar_1.8.1 lattice_0.20-45
## [134] loo_2.5.1 plotrix_3.8-2 fastmap_1.1.0 httr_1.4.4 pkgbuild_1.3.1 glue_1.6.2 shinythemes_1.2.0
## [141] iterators_1.0.14 sass_0.4.4 stringi_1.7.8# END